Ad-aisi  si2 

UNCLASSIFIED 


CLOSING  DEVELOPMENTS  IN  AERODVNAMIC  SIMULATION  HITH 
DISJOINT  PATCHED  MESHES(U)  PEDA  CORP  PALO  ALTO  CA 
C  K  LOMBARD  ET  AL.  28  NOV  8C  AF0SR-TR-87-B782 
F4962B-83-C-8881  F/B  2t/4 


1/1 


NL 


•  o 

UNCLASSIFIED 


AD-A181  512 


WIC  FILE  COFy  @ 


EPORT  DOCUMENTATION  PAGE 

i  t5.  RESTRICTIVE  MARKINGS 


■■■ 


J*  SECURITY  CLASSIFICATION  AUTI^^V 


2b.  OCCLASSIFICATION/DOWNGRA 


4.  PERFORMING  ORGANIZATION  R 


6a.  NAME  of  ferforming  organization 

PEDA  Corporation 


6c.  AOORESS  (City.  Slat*  and  ZIP  Coda) 

4151  Middlefield  Road,  Suite  7 
Palo  Alto,  CA  94303 


|j 


riv 


•a.  NAME  OF  FUNDING/SPONSORING 
ORGANIZATION 

AFOSR 


Ecuador  ESS  (City,  State  and  ZIP  Code) 

Bolling  AFB  DC  20332-6448 


8b.  OFFICE  SYMBOL 
(If  applicable ) 

NM 


3.  OISTRIBUTION/AVAILABILITY  of  report 

Approved  for  public  release;  distribution 
uni imited. 


5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 

.1*-  87-0  782 


7a.  NAME  OF  MONITORING  ORGANIZATION 

Air  Force  Office  of  Scientific  Research 


7b.  AOORESS  (City,  State  and  ZIP  CodeJ  ^  ^  J  Q 

Directorate  of  Mathematical  &  Information 
Sciences,  Bolling  AFB  DC  20332-6448 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 


F49620-85-C-0081 


10  SOURCE  OF  FUNDING  NOS. 


It.  TITLE  (Include  Security  Clarification) 

IN  AERODYNAMIC  SIMULATION  WITH  DISJOINT  PATCHE 


12.  PERSONAL  AUTHOR(S) 

ombard,  E.  Venkatapathy ,  J.  Bardina, 


PROGRAM 
E  LEMENT  NO. 

61102F 

MESHES 


PROJECT 

NO. 


WORK  UNIT 
NO. 


liuniTp 


N.  Nagaraj,  J.  Y.  Yan 


1C  DATE  OF  REPORT  (Yr.,  Mo..  Day) 
86  86  11  28 


R.C^-C.  Luh  & 


15.  PAGE  COUNT 


COSATI  COOES 


18.  SUBJECT  TERMS  (Continue  on  reverge  if  necessary  and  identify  by  block  number) 


Algebraic  Grid  Generation,  Upwind  Methods,  Relaxation, 
Approximate  Factorization,  Euler,  Navier-Stokes,  IjVD  Limiter 


IB.  ABSTRACT  rConlinu*  on  reverse  if  necessary  and  identify  by  block^rrtTrnticrj  * 

In  the  first  three  years  of  a  continuing  progranTthe  research  he^  aimed  at  providing  computational 
tools  and  procedures  as  the  building  blocks  for  a  systenp.  to  permit  efficient  solution  and  high  resolution 
capture  of  flow  structure  in  gasdynamic  problems  of  realistically  complex  geometries.  The  reseafchj  has 
yielded  a  comparatively  simple  algebraic  procedure  for  constructing  two  and  three  dimensional  geometry 
fitted  base  level  composite  meshes  in  quadrilateral  patches.  The  method  provides  complete  control  of 
coordinate  distribution  and  gradient  on  all  patch  boundaries  which  may  include  slope  discontinuities. 
A  robust  upwind  implicit  method  (CSCM)  was  the  basis  to  solve  the  multidimensional  pseudo  time 
dependent  Euler  or  compressible  Navier-Stokes  equations.  Research  into  solution  algorithms  for  that 
upwind  method  has  yielded  a  more  robust  diagonally  dominant  (DDADI)  approximate  factorization 
that  subsequently  led  to  a  family  of  rapidly  convergent  and  data  storage  and  management  efficient 
relaxation  schemes  in  two  and  three  space  dimensions.  JjThose  operationally  explicit  and  unconditionally 


ao.  distribution/availability  of  abstract 
unclasbipieo/unlimited  B  same  as  RPT.  □  OTIC  USERS  □ 


DD  FORM  1473,  83  APR 


21  ./abstract  security  classification 

UNCLASSIFIED 

22h  TELEPHONE  NUMBER 
l Include  Area  Code! 

22c  OFFICE  SYMBOL 

(202)  767-5026 

NM 

EDITION  OF  1  JAN  73  IS  OBSOLETE 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


19.  Abstract  cont. 

stable  upwind  algorithms  have  led  to  a  simple  robust  boundary  procedure  based  on  interpolation  of 
conservative  variable  data  from  adjacent  patches  overlying  interior  patch  boundaries  where  coordinates 
are  discontinuous.  Results  of  tests  with  reflecting  shock  capture  on  overset  adaptively  refined  mesh 
patches  in  a  2-D  supersonic  inlet  problem  show  for  comparable  accuracy  a  savings  of  about  an  order  of 
magnitude  in  mesh  points  relative  to  uniform  mesh  refinement.  A  simple  algorithm  has  been  created  to 
construct  continuous  contour  plots  of  solutions  on  global  domains  .covered  by  such  a  system  of  multiple 
meshes.  The  3/D  symmetric  GaussfSeidel  implicit  method  of  planes  space  marching  relaxation  algorithm 
has  been  implemented  on  a  system  of  composite  patched  meshes  and  applied  in  the  solution  of  a  multi 
rocket  engine  shrouded  exhaust  problem  that  features  large  pockets  of  separated  base  flow.  ~ 

The  methodology  of  the  flux  difference  split  upwind  schemes  has  been  extended  with  new  locally 
one  dimensional  total  variation  diminishing  (TVD)  schemes  that  provide  robust  second  or  third  order 
spatial  accuracy.  One  of  the  latter  schemes  has  been  implemented  in  a  consistent  second  order  method 
in  both  time  and  space  and  applied  in  a  problem  of  unsteady  strong  shock  diffraction  around  a  model 
aeroassisted  orbital  transfer  vehicle  (AOTV).  Finally  in  a  study  of  the  impact  of  mesh  topology  and 
refinement  on  solution  accuracy  we  have  shown  how  local  patches  of  refined  overset  mesh  can  be  used 
to  tie  a  system  of  overlapping  composite  mesh  together  to  give  a  hybrid  topology  with  superior  balance 
in  flow  resolution. 


IIMOT! 


SECURITY  CLASSIFICATION  OF  THIS  PAGI 


AFOSK-m.  8  7-0  782 


Closing  Developments  in  Aerodynamic  Simulation 
with  Disjoint  Patched  Meshes 


Annual  Report 
F49620-85-C-0081 
01  Jun  1985  through  31  Aug  1986 


Charles  K.  Lombard 
Principal  Investigator 
PEDA  Corporation 


Accesion  For 

NTIS  CRA&I 
DTIC  TAB 
Unannounced 
Justification 


By . . 

Oist.  ibutior,  / 


Availability  Codes 

I  AvaiPand/or 
i  Special 


87  6 


2  1  9 


Closing  Developments  in  Aerodynamic  Simulation 
with  Disjoint  Patched  Meshes 

Charles  K.  Lombard,  Ethiraj  Venkatapathy,  Jorge  Bardina, 

N.  Nagaraj,  Jaw  Yen  Yang  and  Raymond  C.-C.  Luh 
PEDA  Corporation 

Joseph  Oliger 
Stanford  University 

1.  Introduction 

Recent  and  projected  developments  in  supercomputers,  numerical  grid  gen¬ 
eration  techniques  and  computational  algorithms  for  the  compressible  Euler  and 
Navier-Stokes  equations  portend  a  major  revolution  in  the  manner,  pace,  cost  of 
design  and  the  resulting  performance  of  aerodynamic  systems.  To  realize  these  po¬ 
tential  benefits,  certain  closing  developments  in  computational  technique  must  be 
made  in  order  to  effect  a  highly  accurate,  reliable,  efficient  and  productive  simula¬ 
tion  environment  for  aerodynamic  design  analysis. 

A  primary  need  of  the  developments  is  to  achieve  the  capability  for  a  user  to 
easily,  rapidly  and  accurately  perform  flowfield  calculations  among  problems  of  dis¬ 
parate  and  realistically  complex  geometries.  The  natural  approach  to  realizing  this 
objective  with  comparatively  straightforward  extensions  of  existing  finite  difference 
computational  technology  is  through  the  use  of  systems  of  quadrilateral  patched 
meshes. 

Such  systems  can  be  either/both  composite  (joined)  or  overset  (disjoint).  In 
the  former  case  adjacent  patches  share  a  common  boundary,  or  at  least  parallel 
boundaries  in  the  case  of  mesh  patch  overlap  for  purposes  of  applying  numerical 
boundary  conditions.  With  composite  meshes,  patch  boundaries  are  piecewise  fitted 
to  segments  of  physical  or  computational  boundaries  or  embedded  flow  structures. 
As  shown  by  Lombard,  et  al1,  composite  mesh  systems,  that  may  have  numerically 
useful  properties  of  geometric  continuity  across  patch  boundaries,  admit  topolog¬ 
ically  singular  global  meshes  that  have  the  capability  to  connect  computational 
regions  of  great  (really  any)  geometric  complexity.  However,  situations  exist  where 
multiple  mesh  topologies,  each  naturally  related  to  some  different  piece  of  geome¬ 
try  or  flow  structure,  offer  greater  flexibility  and  accuracy  than  composite  meshing 
alone.  Examples  involve  multiple  bodies  that  may  have  relative  motion  or  weak 
shocks  that  bear  little  geometric  relationship  to  boundaries  of  the  flow  domain.  In 
such  cases  systems  of  disparately  oriented  and  at  least  partially  overset  grids  as 
proposed  by  Berger  and  Oliger3  allow  arbitrarily  high  resolution  of  all  features  of 
the  flowfield. 

To  make  efficient,  productive  use  of  patched  meshing  strategies  requires  a  body 
of  new  computational  tools  and  methodology  that  are  the  objectives  of  the  present 
research.  The  needed  factors  are:  (1)  a  simple  procedure  for  generating  patched 


computational  meshes  with  freedom  of  point  and  gradient  specification  on  all  patch 
boundaries;  (2)  improved  upwind  algorithm/numerical  boundary  condition  proce¬ 
dures  for  semi-autonomous  implicit  but  unconditionally  stable  conservative  coupling 
of  solutions  on  a  system  of  multiple  patched  meshes;  and  (3)  computer  graphics, 
particularly  a  simple  algorithm  for  constructing  contour  plots  on  systems  of  over¬ 
set  patched  meshes.  The  glue  that  is  to  bind  these  tools  together  in  a  simulation 
requiring  minimum  human  intervention  to  adapt  to  new  configurations  is  a  flexible 
parameter  controlled  multiple  mesh  data  structure.  An  important  objective  of  the 
program  is  to  test  the  evolving  techniques  in  appropriate  problems. 

2.  Research  Accomplishments 

This  annual  report,  the  first  under  contract  F49620-85-C-0081,  summarizes  the 
research  accomplishments  of  the  first  three  years  of  a  continuing  program.  In  the 
first  year  of  the  program  the  emphasis  was  placed  on  the  most  crucial  and  challeng¬ 
ing  objectives  -  patched  grid  generation  and  robust  upwind  algorithm/boundary 
procedures  for  rapid  relaxation  on  multiple  meshes.  The  first  year’s  effort9  sufficed 
to  create  some  needed  tools  of  algebraic  grid  generation  and  to  establish  an  opera¬ 
tionally  explicit  but  unconditionally  stable  upwind  algorithm/numerical  boundary 
condition  procedure  for  systems  of  patched  meshes.  The  techniques  were  tested  in 
comparatively  simple  problems. 

In  the  second  year’s  effort4  the  techniques  were  implemented  in  2-D  and  tested 
against  internal  aerodynamics  problems  of  moderate  geometric  complexity.  The 
upwind  algorithm  was  extended  from  two  to  three  spatial  dimensions.  A  simple, 
effective  algorithm  for  contour  plotting  on  domains  covered  by  multiple  patched 
meshes  was  created  to  exhibit  results. 

In  the  third  year  the  3-D  upwind  implicit  relaxation  algorithm  was  implemented 
on  a  system  of  composite  patched  grids.  The  techniques  in  two  and  three  space 
dimensions  were  tested  against  realistic  internal/external  aerodynamics  problems 
of  challenging  geometric  complexity.  Within  the  framework  of  one  of  the  problems, 
a  continuing  study  of  alternate  composite  and  overset  patched  meshing  topologies 
demonstrated  a  strategy  of  superior  balance  in  mesh  point  utilization  for  overall 
resolution  of  flow  structures  on  the  domain. 

Also  in  the  third  year,  new  implicit  TVD  schemes  were  constructed  to  improve 
the  nonlinear  stability  of  second  and  third  order  accurate  spatial  differencing.  One 
of  the  latter  was  implemented  in  a  consistent  second  order  method  in  both  time  and 
space  and  applied  in  an  unsteady  strong  shock  diffraction  problem  around  a  model 
hypersonic  space  transportation  vehicle  (AOTV). 

Algebraic  Grid  Generation 


The  concept  of  patched  meshing  in  which  complex  domains  are  broken  up  into 
many  geometrically  regular  and  topologically  rectangular  subdomains  leads  natu¬ 
rally  to  the  use  of  efficient  algebraic  techniques  for  the  construction  of  the  individual 


mesh  patches.  To  obtain  the  desired  smoothness  properties  over  the  global  mesh 
in  the  vicinity  of  patch  boundaries,  a  technique  that  permits  specification  of  point 
distribution  and  gradient  on  all  boundaries  was  devised.  The  technique  -  termed 
generalised  transfinite  interpolation6  -  makes  use  of  a  parameterized  general  cubic 
polynomial  for  the  coordinate  curves.  Regularity  of  the  mesh  is  obtained  by  em¬ 
ploying  continuous  distributions  of  the  parameters  of  the  curves  within  judiciously 
chosen  bounds  based  on  analysis.  Stretching  functions  such  as  that  of  Vinokur6  are 
used  to  distribute  points  and  blending  functions  are  used  to  distribute  parameters 
of  the  curves  between  lateral  boundaries. 

A  novel  feature  of  the  technique  is  the  introduction  of  the  corner  singularity 
from  analysis  to  govern  distribution  of  points  and  parameters  in  the  vicinity  of 
boundary  slope  singularities.  At  such  points,  the  method  thus  obtains  the  desired 
properties  of  mesh  smoothness  to  the  interior.  A  global  mesh  solution  obtained  with 
the  method  for  a  backward  step  problem  is  shown  in  Figure  1.  Here  the  solution 
was  generated  in  two  patches,  one  containing  the  exterior  corner  and  the  other  the 
interior  corner.  The  solution  was  matched  analytically  at  the  patch  interface. 

Early  in  the  second  year,  in  the  process  of  attempting  to  apply  the  generalized 
transfinite  interpolation  technique  in  a  variety  of  2-D  problems,  it  became  evident 
the  method  was  too  sensitive  to  parameter  selection  among  too  many  options,  was 
confusing  and  ultimately  required  too  much  artistry  to  meet  the  objectives  of  sim¬ 
plicity  and  user  friendliness  set  for  the  products  of  the  research.  Further,  the  lack 
of  an  analytical  solution  to  corner  problems  blocked  the  straightforward  extention 
of  the  technique  to  3-D. 

With  some  reflection  it  became  evident  to  us  that  the  difficulty  lay  in  trying 
to  accomplish  too  much  in  a  single  step  process.  Rather,  borrowing  the  tools  of 
the  algebraic  technique  and  redefining  the  process  in  multiple  steps  with  interactive 
computer  graphics  sets,  we  could  define  a  straightforward  procedure  to  meet  the 
desired  ends. 

The  approach  that  has  been  implemented  is  in  the  realm  of  two  boundary 
methods,  in  that  one  pair  of  opposite  sides  of  a  patch  is  regarded  as  prescribed  and 
often  includes  a  portion  of  a  physical  boundary.  The  other  pair  of  sides  is  formed 
of  the  left  and  right  limiting  members  of  the  family  of  generalized  cubic  coordinate 
curves  joining  the  initially  given  two  boundaries.  In  either  2-D  or  3-D  the  general 
cubic  coordinate  curve  has  the  simple  form 

1  =  z±+  (n  +  £j<j(u)  +  £2Mu)  (!) 


Equation  (1)  is  a  hermite  interpolation  of  value  (r)  and  gradient  (2)  on  the  two 
boundaries  and  is  parameterized  in  terms  of  u  which  varies  from  zero  to  unity.  The 
scalings  of  o±  and  £2  influence  the  shape  (curvature)  of  the  curve  between  any  pair 
of  end  points.  The  specification  of  a  discrete  set  of  u  values  using  a  generalized 
distribution  function  such  as  that  of  reference  5  defines  the  nodal  intersections  with 
the  other  family  of  coordinate  curves. 

In  our  implementation,  the  left  and  right  limiting  (lateral  bounding)  coordinate 
curves  are  developed  interactively  on  a  graphics  terminal  or  workstation  to  have 
the  desired  configurations.  The  parameters  of  these  lateral  bounding  curves  are 
then  blended  with  polynomial  weighting  functions  to  describe  the  general  cubic 
coordinate  curve  over  the  intervening  also  discretized  interval. 

The  lateral  patch  boundaries  are  essentially  control  devices  that  specify  shape 
and  distribution  to  surrounding  regions.  As  such  they  are  placed  where  needed  - 
at  breaks  in  body  surface  geometry  and  as  terminators  or  transition  guides  from 
regions  of  strong  shape  variation  to  regions  of  very  regular  mesh.  In  fact  the  mesh 
generation  problem,  particularly  for  geometries  with  any  substantial  complexity,  is 
a  problem  of  multiple  length  scales.  The  purpose  of  multiple  patching  is  to  isolate 
regions  of  comparable  scales  and  on  which  subdomains  the  solution  is  comparatively 
regular  and  can  be  conveniently  fit  by  simple  functions. 

Once  a  primary  grid  is  generated  by  the  technique  described  above  it  can  be 
interactively  improved  by  modifying  parameter  blending  and  point  distribution  in¬ 
cluding  point  redistribution  along  the  alternate  family  of  coordinate  lines  implicitly 
defined  by  the  nodes  on  the  cubic  coordinate  curves.  The  latter  operation  is  in 
the  spirit,  if  not  the  detailed  implementation,  of  a  two  step  generalized  transfinite 
interpolation. 

Another  secondary  operation  that  we  employ  is  the  modification  of  coordinate 
lines  in  the  vicinity  of  a  boundary  to  smoothly  enforce  local  normality.  The  latter 
operation  like  all  the  procedures  has  been  programmed  as  a  convenient  tool  requiring 
minimal  input  to  apply  at  a  boundary.  Finally,  a  parameterized  tension  spline  that 
provides  an  analytical  description  of  a  curve  amongst  discrete  data  is  a  tool  that 
has  proved  useful  in  the  latter  operation,  in  effecting  redistribution  of  points  along 
coordinate  curves  of  either  family  and  for  fitting  numerically  specified  boundary 
data.  In  Figures  2  and  3  are  shown  respectively  typical  examples  of  2-D  and  3-D 
mesh  planes  generated  with  the  simplified  technique. 

Upwind  Implicit  Relaxation  Algorithm/Boundary  Procedures 

Under  the  contract  we  have  devised  a  new  single  level  operationally  explicit 
but  effectively  implicit  algorithm  for  gasdynamics.  The  algorithm  is  particularly 
appropriate  for  multiple  patch  mesh  systems  because  each  solution  sweep  operation 
on  any  patch  is  decoupled  from  any  other.  Thus  the  method  is  not  only  very  storage 
efficient  and  simple  to  program  including  the  coupling  at  patch  boundaries  but,  also, 
can  make  excellent  use  of  parallel  computing  in  several  straightforward  ways. 


Previously  the  Beam- Warming  factored  implicit  algorithm7  with  the  Baldwin- 
Lomax  thin  layer  viscous  approximation*  has  provided  the  basis  for  two  similar 
space  marching  (PNS)  procedures9,10  for  the  compressible  Navier-Stokes  equations. 
These  PNS  methods  which  are  highly  efficient  -  requiring  half  the  data  storage  and 
a  small  fraction  of  the  computer  time  of  two  level  time  dependent  methods  -  have 
proven  effective  for  flows11*12  with  favorable  streamwise  pressure  gradient  or  with 
relatively  small  adverse  pressure  gradients.  However,  in  the  presence  of  a  strong  ad¬ 
verse  pressure  gradient  such  as  occurs  in  a  wing  or  fin  root  regions  the  contemporary 
PNS  methods  suffer  numerical  stability  problems  and  may  infer  streamwise  separa¬ 
tion  even  where  separation  doesn’t  occur13.  In  such  unseparated  (perhaps  weakly 
separated)  regions,  numerical  stability  has  been  maintained  at  the  price  of  employ¬ 
ing  large  amounts  of  artificial  viscosity  with  a  resulting  loss  in  predictive  accuracy 
and  knowledge  of  the  actual  state  of  the  flow.  Where  strong  streamwise  separation 
occurs  the  methods  are  unstable  and  cannot  proceed.  Particularly  for  the  increas¬ 
ingly  relevant  laminar  flow  situation  that  will  be  encountered  at  very  high  altitude 
by  aerodynamic  systems  such  as  orbital  transfer  vehicles  (OTV’s)  and  space  shut¬ 
tles,  streamwise  separation  becomes  a  likely  occurrence14  in  compression  corners 
associated  with  canopies,  pods,  flared  bodies,  wing  or  fin  roots  and  deflected  con¬ 
trol  surfaces.  Thus  a  more  general  technique  is  needed  that  is  inherently  stable  for 
all  types  of  upstream  influence.  At  a  minimum  the  mixed  elliptic-hyperbolic  prob¬ 
lem  requires  global  iteration,  preferably  with  type  dependent  differencing.  More 
background  on  this  problem  area  is  given  in  the  first  annual  report3. 

New  Universal  Single  Level  Scheme  CSCM-S 

The  CSCM  flux  difference  eigenvector  split  upwind  implicit  method16,16,17  for 
the  inviscid  terms  of  the  compressible  Navier-Stokes  equations  provides  the  natural 
basis  for  an  unconditionally  stable  space  marching  technique  through  regions  of 
subsonic  and  streamwise  separated  flow.  In  such  regions  the  split  method  can  be 
likened  to  stable  marching  of  each  scalar  characteristic  wave  system  in  the  direction 
of  its  associated  eigenvalue  (simple  wave  velocity).  In  supersonic  flow,  where  all 
eigenvalues  have  the  same  sign,  the  method  automatically  becomes  equivalent  to  the 
referenced  PNS  techniques  based  on  the  Beam- Warming  factored  implicit  method 
with  the  Baldwin-Lomax  thin  layer  viscous  approximation. 

Compared  to  contemporary  central  difference  methods,  the  CSCM  character¬ 
istics  based  upwind  difference  approximation  with  its  inherent  numerical  stability 
leads  to  greatly  reduced  oscillation  and  greater  accuracy  in  the  presence  of  captured 
discontinuities  such  as  shocks,  contacts  and  physical  or  computational  boundaries. 
The  correct  mathematical  domains  of  dependence  that  correspond  with  physical 
directions  of  wave  propagation  are  coupled  with  well  posed  characteristic  boundary 
approximations16  naturally  consistent  with  the  interior  point  scheme.  The  result  is 
a  faster  sorting  out  of  transient  disturbances  and  substantially  more  rapid  conver¬ 
gence  to  the  steady  state.  The  splitting  and  the  associated  time  dependent  implicit 
method  have  been  described  in  detail  in  references  (15)  and  (17)  for  quasi  1-D  and 


2-D  planar  or  axisymmetric  flow. 

In  the  following,  we  will  sketch  the  differences  between  the  time  dependent 
method  and  the  new  space  marching  technique  which  we  designate  CSCM-S.  The 
discussion  will  begin  with  the  quasi  1-D  inviscid  formulation,  present  some  results 
elucidating  the  properties  and  performance  of  the  method,  then  give  additional  de¬ 
tails  entering  into  multidimensional  inviscid  and  thin  layer  viscous  procedures  and 
present  some  2-D  solutions  obtained  with  the  new  single  level  scheme  in  problems 
solved  previously17  with  the  time  dependent  method.  Lastly  \  sketch  a  3-D  im¬ 
plicit  method  of  planes  algorithm  and  give  some  results  for  an  axisymmetric  nozzle 
flow  over  a  backward  facing  step. 

Quasi  1-D  Formulation 

The  general  jth  interior  point  difference  equations  for  the  time  dependent 
CSCM  upwind  implicit  method  for  the  inviscid  advection  terms  is 

(J  +  A+V  +  A- A)%  =  -A+A?)  n  ,  -  A~Aq)n  (2) 

J  -1  j 

where  V  and  A  are  backward  and  forward  spatial  difference  operators.  Here  q  is 
the  conservative  dependent  variable  vector  and  F  is  the  associated  flux  vector.  In 
the  notation  the  interval  averaged  matrices  between  node  points  j  and  j  +  1  are 
indexed  j.  For  simplicity,  the  right  hand  side  of  equation  (2)  is  written  for  the  first 
order  method.  Higher  order  methods  in  space  are  given  with  results  in  references 
15  and  17.  In  equation  (2)  the  CSCM  flux  difference  splitting  is 

(A+  +  A~)Aq  =  A  F+  +  A  F~  =  A  F  (3o) 

with 

A±  =  (AfTI±T~1M~1)A  =  A±A  (3b) 

and 

=  ^(/±sgn(I))  (3c) 

exhibiting  the  similarity  transformation  that  diagonalizes  the  constructed15,17  flux 
difference  Jacobian  A.  Here  A  is  a  diagonal  matrix  of  the  interval  averaged  eigen¬ 
values  that  through  the  truth  function  diagonal  matrices  7*  make  the  decisions 
about  directions  of  characteristic  wave  propagation  and  whether  or  not  to  send 
signal  to  the  solution  point.  Thus  in  equation  (2)  the  piece  of  the  flux  difference 
splitting  A+Ag)y_i  represents  the  convection  of  characteristic  wave  contributions 
in  the  positive  coordinate  direction  from  grid  point  j  —  1  to  solution  point  j  and  A~ , 
in  the  negative  direction  from  j  +  1  to  j.  As  the  result  of  incorporating  multiplica- 
tively  the  (local)  time  step  (for  pseudo  time  relaxation)  and  the  spatial  (divided) 
differences  in  the  matrices,  the  numerical  eigenvalues  are  Courant  numbers  for  the 
characteristic  waves  whose  speeds  are  u,  u  4-  c  and  u  —  c,  with  c  the  sound  speed. 


Central  to  its  accurate  shock  capturing  capability,  the  CSCM  conservative  flux 
difference  splitting  has  the  Roe18  “property  U”  embodied  in  equation  3a. 

With  6q  =  qn+1  —  qn,  equation  (2)  defines  a  two  level  linearized  coupled  block 
matrix  implicit  scheme  that  can.  be  solved  by  a  block  tridiagonal  procedure.  In 
reference  (17)  a  new  (DDADI)  approximately  factored  alternating  sweep  bidiagonal 
solution  procedure  for  equation  (2)  is  presented  that  is  shown  to  be  very  robust 
and  is  operationally  explicit,  i.e.  requires  only  a  decoupled  sequence  of  local  block 
matrix  inversions  rather  than  the  solution  of  the  coupled  set.  For  the  forward  sweep 
the  bidiagonal  solution  procedure  can  be  written 

(I  +  A+-  A~)6q*j-  =  RHS  +  A+6q*3._ x  (4) 

Fc  the  linear  problem,  i.e.  constant  coefficient  case  of  stability  analysis,  equation 
(4)  is  equivalent  to  the  single  level  space  marching  procedure 

(/  +  A+  -  A-)6q%-  =  A+q*j-i  ~  "  A~Aq)n  (5) 

J  J 

Nonlinearity  enters  in  the  single  level  space  marching  form  (5)  in  that  at  each  step  of 
the  forward  sweep  the  matrices  A+  are  averaged  between  0*y_i  and  g”  rather  than 
homogeneously  at  the  old  iteration  level  n.  Similarly,  a  companion  backward  space 
marching  sweep  that  is  symmetric  to  equation  (5)  and  that  is  intimately  related  to 
the  backward  sweep  of  the  alternating  bidiagonal  algorithm  of  reference  (17)  is 

U  +  A+-  2-)%  =  +  A-,;  -  A~,”  +  J  (6) 

The  method  given  by  equations  (5)  and  (6)  is  von  Neumann  unconditionally  star 
ble  for  the  scalar  wave  equation.  The  analysis  shows  the  significance  of  DDADI 
approximate  factorization  in  rendering  both  the  forward  and  backward  sweeps  sep¬ 
arately  stable  regardless  of  eigenvalue  sign.  Consequently  as  the  local  Courant 
number  becomes  very  large,  the  robust  method  becomes  a  very  effective  (symmetric 
Gauss-Seidel)  relaxation  scheme  for  the  steady  equations,  a  fact  which  substantially 
contributes  to  the  very  fast  performance  that  will  be  demonstrated. 

At  a  right  computational  boundary  on  the  forward  sweep  we  solve  the  charac¬ 
teristic  boundary  point  approximation17 

q#+1  =  q  #  and  at  a  left,  on  the  backward  sweep 


(7) 


Following  the  solution  of  equations  (7)  and  (8)  the  conservative  state  vector  is 
iteratively  corrected16  to  maintain  the  accuracy  of  prescribed  boundary  conditions 
while  not  disrupting  the  representation  of  the  computed  characteristic  variables 
running  to  the  boundary  from  the  interior.  Analysis  of  a  model  system  with  upwind 
differenced  scalar  equations  and  coupled  boundary  conditions  was  related  to  the 
linearized  bidiagonal  scheme17  by  Oliger  and  Lombard10;  the  analysis  also  strongly 
supports  the  numerically  confirmed  robust  stability  of  the  present  nonlinear  method 
for  gasdynamics.  A  useful  result  of  reference  19  that  simplifies  the  procedure  of 
reference  17  is  that  on  the  forward  sweep  there  is  no  need  for  a  predictor  step  at 
the  left  boundary  J  =  1;  thus,  the  solution  sweep  begins  at  J  =  2.  Similarly,  the 
backward  sweep  begins  at  J  =  N  —  1. 

With  the  updating  at  each  step,  where  in  equation  (6)  6qj  =  ?y+1  —  it  is 
clear  that  the  symmetric  pair  of  equations  (5)  and  (6)  serve  to  advance  the  solution 
two  pseudo  time  (iteration)  levels;  whereas,  the  linear  alternating  bidiagonal  sweep 
algorithm  of  reference  (17)  advances  the  solution  only  one  level.  To  maintain  con¬ 
servation  to  a  very  high  degree,  in  single  sweep  marching  in  supersonic  zones  we 
iterate  (at  least)  once  locally  at  each  space  marching  step.  The  local  iteration  serves 
to  make  the  eigenvectors  in  the  coefficient  matrices  consistent  with  the  advanced 
state  and  thus  provides  improved  accuracy  for  the  nonlinear  system.  It  appears 
effective  to  do  this  inner  iteration  everywhere,  i.e.  in  both  subsonic  and  supersonic 
regions,  as  the  number  of  global  iteration  steps  to  convergence  with  two  inner  iter¬ 
ations  has  been  found  reduced  by  a  factor  of  three  to  four.  Since  the  computational 
work  per  two  steps  is  about  the  same  for  the  single  level  and  two  level  schemes,  and 
beyond  the  fact  that  one  saves  a  level  of  storage  in  the  space  marching  algorithm, 
the  question  arises:  Can  one  get  solutions  in  less  computational  work  through  faster 
convergence  with  the  nonlinear  space  marching  algorithm? 

One  Dimensional  Results 

First,  we  present  results  for  supersonic  flow  with  no  shock  in  Shubin’s  diverging 
nozzle.  In  purely  supersonic  zones,  the  experience  with  the  present  method  is  that 
the  solution  can  be  marched  accurately  in  one  global  iteration,  as  ought  to  be  the 
case.  Figure  4  shows  the  exact  solution  (in  solid  line)  and  the  computed  result 
from  the  first  forward  sweep.  It  is  evident  that  the  method  correctly  predicts 
the  solution  to  plotting  accuracy  in  one  forward  sweep.  With  subsequent  sweeps 
the  error  (the  difference  between  the  exact  and  the  computed  solution)  reduces  to 
machine  accuracy  in  less  than  three  global  iterations.  In  fact,  by  increasing  (from 
two)  the  number  of  inner  iterations  on  the  solution  procedure  at  each  space  marching 
step,  convergence  to  prescribed  accuracy  can  be  guaranteed  in  one  forward  sweep. 
This  is  also  true  of  contemporary  locally  linearized  unsplit  methods  in  supersonic 
flow. 

With  the  globally  iterative  nonlinear  space  marching  formulation,  early  expe¬ 
rience  in  two  quasi  1-D  nozzle  problems  with  mixed  supersonic-subsonic  zones  is 
that  solutions  are  obtained  in  roughly  an  order  of  magnitude  fewer  iteration  steps 


than  had  been  required  with  the  previously  fast  pseudo  time  dependent  technique 
and  block  tridiagonal  solving. 

The  two  nozzle  problems  which  are  described  and  solved  by  Yee,  Beam  and 
Warming20  and  solved  with  the  CSCM  time  dependent  technique  in  references  (15) 
and  (16)  are  Shubin’s  diverging  nozzle  flow  and  Blottner’s  converging-diverging 
nozzle  flow.  Both  problems  involve  unmatched  overpressures  at  the  outflow  which 
result  in  internal  shock  terminated  supersonic  zones  and  subsonic  outflow.  For 
the  experiments  involving  flow  of  mixed  type  the  same  initial  data  given  by  Yee, 
Beam  and  Warming  -  a  linear  interpolation  between  inflow  and  outflow  values  for 
effectively  exact  solutions  of  the  problems  -  is  used  that  was  used  previously  with 
the  time  dependent  approach. 

For  flows  of  mixed  type,  in  Figures  5  and  6  respectively,  results  are  shown  for 
successive  forward  and  backward  sweeps  for  five  global  iteration  steps  with  Shubin’s 
and  Blottner’s  nozzle  flows.  In  both  cases,  the  exact  solution  as  given  by  Yee,  Beam 
and  Warming  is  shown  in  solid  line  and  the  present  computational  results  solved 
on  a  51  point  mesh,  in  boxes.  Blottner’s  nozzle  flow  is  shown  converged  after  10 
global  iteration  steps.  There  is  substantial  evidence  in  other  results  not  shown  that 
with  further  work  the  number  of  global  iterations  required  to  compute  flows  such 
as  Blottner’s  can  be  reduced  by  a  factor  of  two,  to  about  five. 

In  Figure  7,  we  show  a  subcritical,  i.e.  completely  subsonic,  flow  solution 
computed  in  only  two  global  iteration  steps  for  the  Blottner  nozzle  geometry  with 
different  inflow  conditions.  Here  the  exact  analytical  solution  derived  by  Venkata- 
pathy  is  shown  in  solid  line  and  our  computed  results  in  boxes. 

The  alternating  direction  sweeps  in  our  method  have  been  derived  directly  out 
of  theory  for  solving  the  implicit  set  of  difference  equations.  However,  one  can  see 
mechanistically  (numerically  speaking)  that  omitting  the  backward  sweep  from  the 
pair  and  globally  iterating  only  with  the  forward  sweep  equation  (5)  will  result  in 
permitting  the  influence  of  a  subsonic  outflow  boundary  (or  interior  disturbance)  to 
propagate  upstream  only  one  grid  point  per  global  iteration.  In  such  a  case,  which 
relates  to  other  global  iteration  methods  found  in  the  literature  and  that  also  sweep 
only  in  the  main  flow  direction,  the  rate  of  convergence  is  greatly  inhibited  relative 
to  symmetric  sweeping  by  a  factor  of  order  roughly  the  number  of  grid  points  in  the 
subsonic  zone.  Mathematically,  this  inhibition  is  the  result  of  the  failure  to  include 
in  the  implicit  process  the  effect  of  the  eigenvectors  governing  upstream  influence 
but  instead  treating  these  waves  explicitly  with  effective  CFL  unity. 

In  Figure  8  we  illustrate  the  progress  of  the  transient  solution  to  the  subcritical 
nozzle  problem  after  15  forward  sweeps,  with  the  backward  sweeps  omitted.  One 
can  clearly  see  that  the  wave  influence  of  the  outflow  boundary  has  progressed  only 
15  mesh  points  forward  of  the  outflow  boundary.  In  Figure  9  the  transient  solution 
is  shown  after  60  steps  which  is  beyond  one  characteristic  transit  time  (equivalent 
to  50  mesh  intervals)  for  the  upwind  wave  to  reach  the  inflow  boundary.  In  Figure 
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10  we  show  the  history  of  the  RMS  error  in  the  primitive  variables.  The  solution  is 
found  to  converge  to  roughly  the  same  RMS  error  after  three  characteristic  times 
(150  steps)  as  the  solution  obtained  with  the  symmetric  alternating  sweep  sequence 
after  only  3  global  iterations. 

Blottner’s  supercritical  nozzle  problem  which  involves  subsonic  inflow  acceler¬ 
ating  through  a  sonic  point  to  a  supersonic  zone  terminated  by  a  shock  to  subsonic 
outflow  is  the  most  computationally  demanding  of  the  test  cases  and  indicates  the 
capability  for  the  method  to  compute  simply  and  consistently  over  the  subsonic 
forebody  and  base  regions  of  blunt  bodies  in  supersonic  flow.  Thus  the  need  for 
separate  time  dependent  codes  will  be  obviated  by  this  new  method. 

Finally,  in  Figures  11  and  12,  we  present  the  convergence  history  for  the  present 
nonlinear  scheme  and  the  linearized  time  dependent  scheme  for  completely  subsonic 
and  supersonic  nozzle  flows.  The  x-axis  shows  the  number  of  iterations  each  scheme 
requires  to  reduce  the  exact  error  to  five  orders  of  magnitude  for  various  Courant 
numbers.  It  is  evident  that  the  present  scheme  converges  extremely  fast  at  all  CFL 
numbers  compared  with  the  method  based  on  the  linearized  block  tridiagonal  solver. 
Two  Dimensional  Formulation 

For  two  dimensional  flow,  assuming  a  marching  coordinate  £,  inviscid  terms 

B+Vn+B~ A„  (9a) 

and 

-B+A^k-x  -B-&r,q)k  (96) 

are  added  to  the  left  and  right  hand  sides  respectively  of  both  the  forward  and  back¬ 
ward  sweep  equations  (5)  and  (6).  For  viscous  flow,  second  centrally  differenced, 
thin  layer  viscous  terms  are  also  added  in  the  rj  direction  as  is  conventionally  prac¬ 
ticed,  e.g.  Steger21.  With  the  terms  for  the  r?  cross  marching  coordinate  direction, 
the  technique  now  becomes  an  implicit  method  of  lines.  Along  each  rj  coordinate 
line,  one  can  solve  the  equations  coupled  with  a  block  tridiagonal  procedure.  Alter¬ 
natively,  a  further  DDADI  bidiagonal  approximate  factorization  can  be  employed  in 
the  rj  direction  and  solved  either  linearly  as  in  reference  (17)  or  nonlinearly  as  here 
in  the  £  direction.  As  shown  in  the  quasi  2-D  numerical  experiments  of  reference 
(17),  DDADI  bidiagonal  approximate  factorization  is  stable  for  viscous  as  well  as 
inviscid  terms.  Finally  in  reference  (17)  there  is  a  relevant  discussion  of  the  reduced 
approximate  factorization  error  that  attends  using  DDADI  in  one  or  more  space 
directions. 

Two  Dimensional  Results 

We  present  results  for  a  45°  -  15°  axisymmetric  transonic  nozzle  flow  previously 
studied  experimentally  by  Cuffel,  Back  and  Massier22  and  computationally  by  Cline, 
Prozan,  Serra  and  Shelton  (all  referenced  in  (22))  and  ourselves17.  In  Figure  13  we 
show  results  after  10  steps  of  an  early  computation  run  at  a  local  CFL  number  of 
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14  with  the  present  first  order  single  level  scheme.  Except  for  the  addition  of  an 
error  correction  procedure1*'  to  counter  numerical  inflow  boundary  condition  drift, 
a  factor  which  has  improved  the  present  solution  in  the  vicinity  of  the  axis,  the 
effectively  converged  results  found  here  are  the  same  as  those  given  for  the  two 
level  scheme  in  reference  17.  (As  long  as  the  problem  has  a  unique  solution,  the  two 
schemes  must  give  equivalent  results  since  the  right  hand  side  difference  equation 
sets,  including  boundary  approximations,  are  the  same.) 


For  the  solution  given  in  Figure  13,  we  noted  a  very  rapid  rate  of  reduction 
in  residual,  three  orders  of  magnitude  in  ten  steps.  This  compares  with  60  steps 
given  in  reference  (17)  for  the  solution  obtained  with  the  two  level  scheme.  The 
rapid  convergence  found  in  this  transonic  problem  for  the  CSCM-S  method  with 
viscous  terms  provides  the  reasonable  expectation  of  similar  fast  results  to  be  ob¬ 
tained  without  viscous  effects.  Thus  the  method  in  multidimensions  appears  to  have 
attractive  potential  for  an  improved  transonic  Euler  solver  as  well  as  Navier-Stokes 
solver. 


Next,  we  present  first  order  inviscid  and  viscous  results  for  an  inlet  problem 
shown  in  Figure  14.  The  pressure  contours  for  the  first  order  inviscid  solutions 
are  shown  in  Figure  15.  Figure  16  shows  the  first  order  viscous  results.  The  vis¬ 
cous  computation  shows  the  presence  of  the  leading  edge  shock.  The  flow  structure 
compares  very  well  with  the  theoretical  (for  the  inviscid  case)  and  other  computa¬ 
tional  results.  In  Figure  17,  the  inviscid  and  viscous  wall  pressure  are  compared 
with  the  exact  solution  (inviscid).  Figure  18  shows  the  convergence  history  of  the 
RMS  residual  of  all  the  conservative  flow  variables  for  the  inviscid  problem  solved 
at  CFL  =  100  with  4  inner  iterations  at  every  axial  location.  For  the  inviscid  case, 
only  forward  marching  was  carried  out  and  backward  marching  was  omitted.  The 
solution  has  converged  for  practical  purposes  at  the  end  of  the  first  sweep.  The 
residual  reaches  machine  accuracy  in  i.0  iterations.  In  reference  23,  we  show  the 
residual  reduction  versus  inner  iteration  number  in  single  sweep  solutions  for  super¬ 
sonic  flow  and  compare  results  with  contemporary  PNS  procedures.  Finally  in  work 
described  in  reference  24  we  have  applied  the  single  level  scheme  to  the  solution  of 
the  coupled  forebody  captured  shock  layer  and  massively  separated  wake  flow  of  a 
hypersonic  axisymmetric  AOTV. 


Three  Dimensional  Method  of  Planes  Algorithm 


In  reference  23  we  presented  a  symbolic  algebra  for  DDADI  approximate  fac¬ 
torization  and  derived  single  level  relaxation  schemes.  The  algebra  is  based  on  the 
implicit  difference  stencil  of  the  implicit  method.  Here  we  will  show  how  the  ap¬ 
proach  can  be  used  to  derive  an  symmetric  Gauss-Seidel  implicit  method  of  planes 
relaxation  algorithm. 


The  unfactored  three  dimensional  linearized  implicit  method  can  be  represented 
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by  the  symbolic  matrix  expression 


B~ 

I  c- 

-A+  - -  D  -  A~ 

-c+  | 

-B+ 


On  the  block  diagonal  the  matrix  D  =  I  +  A+  —  A~  +  B+  —  B~  +  C+  —  C~ .  A 
once  DDADI  approximate  factorization  in  the  £  coordinate  direction  leads  to  the 
expression 


-A+ 


B- 


-B+ 


D"1 


By  analogy  with  the  derivation  of  the  single  level  scheme  for  quasi  1-D  flow  from 
DDADI  bidiagonal  approximate  factorization  we  identify  the  above  expression  with 
the  alternate  space  marching  implicit  method  of  planes  algorithm 

Forward  Sweep 


B~ 

-C+  D  C~ 
-B+ 


Sii  =  RHS\  1 

Backward  Sweep 


B~ 

-C+  D  C~ 
-B+ 


6q»+l  =  rhs[  9;+1 ,  ] 


n+l  „n+2 


In  the  planes  the  coupled  block  matrix  problem  can  be  further  simplified  by  the 
approximate  factorization 


[  -C+  ,  D  ,  C~  }  D-1  [  -B+  D  B~  }6q  =  RHS  (lOo) 

which  leads  directly  to  the  block  tridiagonal  solution  sequence 

[  -C+  D  C~  )Sq *  =  RHS  (10 b) 

(  -B+  D  B~  )6q  =  D6q *  (10c) 
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Three  Dimensional  Results 

This  3-D  space  marching  algorithm  has  been  tested  against  an  axisymmetric 
viscous  flow  problem  of  a  supersonic  RL-10  nozzle  exhausting  over  a  backward  facing 
step  into  a  cylindrical  shroud. 

In  Figure  19  we  show  a  3-D  perspective  view  of  the  wall  surface  mesh  in  the 
quarter  sector  for  which  we  solved.  A  section  of  computational  mesh  in  a  longitu¬ 
dinal  plane  through  the  axis  is  shown  in  Figure  20.  For  such  planes,  Figures  21  and 
22  show  pressure  contour  and  velocity  vector  plots  that  exhibit  the  expected  flow 
structures  -  a  weak  shock  in  the  nozzle  diffuser,  a  strong  expansion  at  the  nozzle 
exit  with  subsequent  recompression  shock  off  the  shroud  wall  and,  lastly,  a  substan¬ 
tial  region  of  streamwise  separated  flow  under  the  backward  facing  step.  The  Mach 
contour  plot  of  Figure  23  for  a  cross  flow  data  plane  in  the  shroud  shows  the  ex¬ 
cellence  with  which  the  method  and,  in  particular  the  lateral  boundary  conditions, 
numerically  maintain  the  axisymmetry  of  the  solution. 

Patched  Mesh  Boundary  Procedures 

In  previous  work  Lombard,  et  al16,17  and  Oliger  and  Lombard10  gave  stable 
implicit  procedures  for  computing  the  solution  at  external  boundaries  of  a  com¬ 
putational  domain.  Those  procedures  generalized  the  work  of  Kenzer  to  matrix 
coupled  linearized  boundary  conditions  complementing  a  set  of  advective  difference 
equations  (associated  with  well  posed  characteristics)  to  the  boundary. 

Under  the  present  contract  we  have  explored  the  problem  of  implicitly  cou¬ 
pling  at  interior  patch  boundaries  the  global  solution  on  a  system  of  multiple  patch 
meshes.  The  approach  we  have  taken  in  this  research  is  numerical  experimenta¬ 
tion  among  a  number  of  boundary  treatment  approximations  to  the  solution  of  the 
continuous  domain  problem.  For  comparison  purposes  with  previous  single  mesh 
results,  numerical  experiments  were  performed  with  the  well  tested  two  level  pseudo 
time  relaxation  CSCM  scheme  at  the  interior  points  of  the  mesh  patches. 

Quasi  2-D  Studies 

Our  first  experiments,  to  be  described  here,  dealt  with  breaking  a  single  co¬ 
ordinate  line  into  segments  and  solving  sequentially  on  each  the  equations  for  a 
quasi  2-D  viscous  compressible  flow.  The  model  problem,  with  which  we  experi¬ 
mented  in  reference  17  for  an  uninterrupted  mesh,  is  a  transient  pipe  flow  resulting 
from  an  initially  nonequilibrated  pressure  between  the  axis  and  wall  boundaries. 
Three  kinds  of  cases  were  run  with  the  two-level  linearized  implicit  procedure;  all 
featured  sequential  solving  on  patches  with  at  least  one  point  of  overlap.  Case  1 
had  frozen  boundary  data,  obtained  from  the  solution  on  neighboring  meshes  in  an 
operationally  explicit  manner.  Specifically  at  left  and  right  first  computed  interior 
points  of  a  patch,  we  solved  the  bidiagonal  equations  respectively 


and 


(/-  A-  +  A+V)«,w_,  =  A+9jv"  2  -  (A+  -  A-)?^"  j  -  A",  "  (116) 

Here  the  symbols  n  indicate  the  “frozen”  boundary  data  coming  from  the  solution 
on  the  interior  of  an  adjoining  and  partially  overlying  patch  which  was  at  the  old 
iteration  level  in  the  global  procedure.  In  Case  1  the  solution  on  all  the  grids  was 
effectively  updated  at  the  same  time. 

Case  2  featured  reverse  sequential  cycling  of  the  two  level  time  dependent  solu¬ 
tion  procedure  through  the  grids  on  alternate  global  steps.  In  the  forward  sequence 
the  right  interior  patch  boundaries  were  (a)  frozen  or  (b)  computed  with  one-sided 
characteristic  boundary  conditions  obviating  any  change  in  the  characteristic  data 
from  outside  (to  the  right  of)  a  patch.  The  left  interior  patch  boundaries  inherited 
implicit  data  (6?)  from  the  solution  on  the  computed  patch  to  the  left.  In  the 
backward  sequence  all  the  roles  were  reversed  including  the  directions  of  forward 
elimination  and  back  substitution  in  the  tridiagonal  matrix  inversion  procedure. 

Case  3  featured  cycling  through  the  patches  in  the  predictor  (forward  elimi¬ 
nation)  step  of  the  solution  procedure,  inheriting  implicit  left  boundary  data  as  in 
Case  2.  Then  the  patches  were  cycled  through  in  reverse  order  on  the  corrector 
(back  substitution)  step.  The  result  (save  for  interpolation  if  data  points  of  the 
grids  were  interlaced)  is  identically  equivalent  to  solving  on  an  uninterrupted  single 
patch  mesh. 

The  results  of  the  three  sets  of  experiments  were  comparable  for  local  time 
steps  based  on  constant  CFL  number  up  to  about  50.  Beyond  that  the  rate  of 
convergence  of  the  solutions  for  Cases  1  and  2  diminished  and  at  sufficiently  high 
CFL  number  failed  to  converge.  The  effort  involved  in  Case  2  with  the  two  level 
scheme  was  not  rewarded  relative  to  the  simple  procedure  of  Case  1.  Within  the 
framework  of  the  single-level  symmetric  Gauss-Seidel  implicit  relaxation  scheme, 
which  is  linearly  equivalent  to  the  bidiagonal  scheme  employed  at  the  boundary, 
Case  2a  is  operationally  no  more  difficult  than  Case  1  and  more  closely  approximates 
Case  3.  Case  2b  has  a  consistency  problem  that  inhibits  firm  convergence. 

Figure  24  compares  RMS  residual  (density)  convergence  history  for  Cases  1 
and  3  with  three  mesh  segments  with  two  cells  of  overlap  and  run  at  CFL  25.  As 
one  might  expect,  the  effectively  uninterrupted  mesh  procedure  is  found  to  converge 
(two  to  three  times)  faster  for  about  the  first  50  steps  through  the  major  transient. 
After  that  the  performance  is  comparable.  The  performance  of  the  frozen  boundary 
treatment  Case  1  with  five  segments  at  CFL  25  was  found  to  be  not  materially  worse 
than  with  three  segments,  but  the  performance  degradation  with  increasing  CFL 
was  found  to  proceed  faster  with  increasing  number  of  patches. 

From  the  results  of  the  quasi  2-D  experiments  and  extrapolation  from  experi¬ 
ence  with  the  single  level  scheme,  we  observe  that  the  simple  boundary  procedure 
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with  frozen  conservative  variable  data  taken  from  the  solutions  at  either  iteration 
level  n  or  n  + 1  on  adjacent  meshes  and  coupled  implicitly  by  alternating  direction 
sequential  solving  through  the  patches  has  robust  stability  to  sufficiently  high  CFL 
number  to  yield  a  rate  of  convergence  meeting  our  needs. 

Two  Dimensional  Patch  Boundary  Treatment 

In  earlier  work  solving  the  Euler  equations  on  multiple  patch  meshes  Benek, 
et  al.26  employed  linear  interpolation  of  the  conservative  variables  from  the  interior 
solution  of  one  mesh  to  give  Dirichlet  boundary  data  for  the  other.  With  sev¬ 
eral  points  of  mesh  overlap  at  the  mesh  boundaries,  transonic  solutions  obtained 
with  central  differencing  exhibited  considerable  oscillation  in  the  vicinity  of  shocks 
propagating  through  the  boundary  region.  Eberhardt27  with  the  code  of  Benek, 
et  al.26  attempted  to  reduce  the  oscillation  and  attendant  stability  problems  en¬ 
countered  in  the  vicinity  of  a  bluff  body  shock  intersecting  an  embedded  patch 
boundary  by  introducing  a  characteristic  computed  boundary  point  approximation 
with  scalar  upwind  difference  equations  in  Riemann  variables.  In  his  procedure, 
interpolation  was  performed  only  on  the  variables  whose  characteristics  ran  to  the 
patch  boundary  from  outside  the  computational  domain.  Eberhardt  based  the  de¬ 
cision  about  domain  of  dependence  of  the  characteristics  on  eigenvalues  computed 
within  the  patch  domain.  When  this  decision  was  compatible  with  the  flow,  then 
the  characteristics  boundary  procedure  gave  markedly  superior  results  compared 
to  interpolation  of  conservative  variable  data,  which  leads  to  solution  overspecifica¬ 
tion  in  subsonic  zones.  In  other  cases  where  incorrect  domain  of  dependence  was 
inferred,  the  characteristic  boundary  procedure  was  unstable. 

Here,  based  on  our  quasi  2-D  studies  described  above,  we  give  a  simple,  ro¬ 
bustly  stable  implicit  approach  to  computing  solutions  of  the  conservative  equations 
of  gasdynamic8  on  either  composite  or  overset  meshes.  Without  requiring  special 
flux  conservative  operators,  but  rather,  interpolating  conservative  variable  data  at 
mesh  boundaries,  the  implicit  upwind  method  is  accurate  and  relatively  free  of 
oscillation  where  shocks  intersect  interior  patch  boundaries.  In  a  supersonic  inlet 
problem  with  an  expansion  and  reflected  shocks,  we  exhibit  the  capability  to  conve¬ 
niently  carry  out  with  rapidly  convergent  implicit  methods  for  systems  of  equations 
the  adaptive  refined  meshing  strategy  in  overset  patches  proposed  by  Berger  and 
Oliger2.  Further,  the  test  case  shows  concretely  in  a  realistic  aerodynamic  problem 
the  savings  in  mesh  points  (about  an  order  of  magnitude  here  in  two  dimensions) 
for  similar  accuracy  that  flow  structure  aligned  adaptive  patched  meshing  affords 
compared  to  uniform  grid  refinement. 

The  factors  in  our  approach  are  supported  by  previous  research  by  ourselves  and 
generically  by  others  cited  in  reference  28  and  are  proven  in  numerical  experiments 
reported  here  and  elsewhere22.  We  employ  an  implicit  conservative  upwind  scheme 
CSCM17  with  which  in  the  present  work  we  can  solve  to  either  first  or  second 
order  spatial  accuracy  the  Euler  or  compressible  Navier-Stokes  equations  in  two- 
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dimensional  planar  or  axisymmetric  flows.  The  flux  difference  split  upwind  schemes 
of  generalized  Roe  form  such  as  CSCM  have  a  number  of  qualities  that  make  them 
ideal  for  the  purpose  of  solving  on  discontinuous  patched  mesh  systems. 

First,  conservative  schemes  in  the  Roe  form  satisfy  Roe’s  property  U  that 
guarantees  the  correct  speed  for  captured  shocks.  The  CSCM  scheme  has  been 
tested  in  a  wide  variety  of  internal  and  external  transonic  and  hypersonic  flows17,23 
and  has  been  found  to  capture  strong  and  weak  shocks  accurately  in  location  and 
with  little  oscillation.  The  shock  transition  is  particularly  sharp,  about  two  mesh 
cells  wide,  on  an  aligned  grid;  and  this  factor  will  be  accommodated  as  much  as 
possible  in  our  adaptive  patched  mesh  strategy. 

Second,  in  the  Roe  form,  the  difference  operators  on  conservative  variable  data 
represent  the  effects  of  differencing  to  characteristic  data  only  for  disturbances  prop¬ 
agating  toward  the  given  node  and  reject  the  mathematically  unstable  contribution 
from  disturbances  that  may  be  propagating  downwind  of  the  node.  The  one  sided 
upwind  difference  operation  represents  identically30  the  (split)  conservative  partial 
flux  difference  between  the  nodes.  To  the  extent  that  the  data  from  an  adjacent 
mesh  is  consistent  with  the  solution,  the  associated  upwind  partial  flux  difference 
to  that  data  will  serve  to  provide  at  convergence  consistency  of  the  partial  flux 
convective  into  the  given  mesh  from  its  neighbor,  and  vice  versa  for  signals  of  the 
other  eigenvalue  sign.  Thus  the  method  acts  within  truncation  error  to  provide  the 
similar  continuity  of  the  flux  tensor  among  patched  grids  that  exists  in  the  physical 
solution  across  shocks  and  that  would  obtain  on  a  single  grid  alone.  The  correct 
domain  of  dependence  coupled  with  the  well  posed  characteristic  boundary  point 
approximations17  tend  not  to  support  oscillatory  disturbances  but  convect  them 
out  of  the  flow  domain. 

Third,  one  sided  difference  interval  averaged  eigenvalues  let  majority  rule  de¬ 
termine  the  direction  of  local  signal  propagation.  When  applied,  as  we  do,  to  a 
difference  operator  between  boundary  data  obtained  from  an  adjacent  mesh  and 
the  local  mesh  solution  point,  the  data  of  both  meshes  participate  in  making  the 
decision  as  to  whether  an  incoming  signal  is  being  sent.  Both  in  concept  and  from 
our  experience,  this  factor  seems  to  overcome  the  inter-mesh  communication  diffi¬ 
culty  experienced  by  Eberhardt27  with  his  characteristic  boundary  procedure. 

Lastly,  with  the  CSCM  difference  equations,  with  diagonally  dominant  approx¬ 
imate  factorization17,23  that  retains  on  the  diagonal  the  contributions  from  both 
sets  of  eigenvalues  in  what  is  effectively  an  absolute  value  Jacobian  matrix,  we  can 
solve  the  equations  either  with  two  data  level  linearized  block  implicit  methods17 
or  with  a  single  data  level  relaxation  technique23  that  is  substantially  more  rapidly 
convergent  than  the  linearized  implicit  procedures.  As  can  be  inferred  from  the 
theory  and  numerical  experiments  of  reference  23,  the  use  of  DDADI  on  the  solu¬ 
tion  point  while  differencing  effectively  explicitly  to  data  obtained  from  an  adjoining 
mesh  (which  may  be  at  either  the  old  (n)  or  new  (n+1)  level)  is  unconditionally  sta¬ 
ble.  When  the  solutions  on  the  patched  meshes  are  alternately  updated  using  either 
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the  linearized  implicit  or  relaxation  methods  the  global  procedure  is  implicit.  As 
will  be  shown  here,  the  robust  stability  of  the  global  procedure  has  been  confirmed 
to  approximately  coarse  mesh  CFL  100. 

One  point  that  has  not  been  touched  on  is  the  form  of  interpolation  that  we 
use.  Differencing  to  interpolated  data  is  equivalent  to  a  weighted  sum  of  differences 
operating  on  the  interpolants.  It  is  intuitive  that  for  robust  stability  each  of  these 
assumed  upwind  differences  ought  to  be  well  posed.  This  implies  that  the  inter¬ 
polation  weights  should  all  be  positive  and  the  domains  of  dependence  of  all  the 
interpolants  should  be  outside  (i.e.  on  the  assumed  side)  of  the  solution  point  with 
respect  to  its  mesh  interior.  Neither  of  these  properties  was  shared  by  the  data 
interpolation  schemes  used  by  Benek,  et  al.20  or  Eberhardt27. 

Finally,  extending  a  direction  of  Benek,  et  al.20,  we  do  not  compute  on  sections 
of  coarser  meshes  underlying  patches  of  overset  refinement.  Our  data  structure  and 
automated  procedures  for  the  consequent  partitioning  of  meshes  are  described  in 
reference  31.  The  partitioning  concept  in  which  coordinate  lines  of  a  patch  have 
differing  (index)  lengths  in  computational  space  also  leads  to  useful  possibilities 
including,  as  will  be  shown  here,  fitting  mesh  patches  obliquely  to  boundaries,  e.g. 
to  sharply  capture  reflecting  shocks. 

Interior  Boundary  Treatment  -  First  Computed  Point  Formulation 

In  Equations  (2)  -  (8)  we  sketched  the  development  of  the  single  level  relaxation 
algorithm  and  the  treatment  at  external  boundaries  of  the  global  computational 
domain.  Here  we  describe  the  treatment  of  an  interior  patch  boundary.  The  left 
hand  side  of  equation  (2)  has  the  tridiagonal  structure 

A+  ,  I+A+-A-  ,  A~  (12) 


In  relation  (12),  the  central  block  which  we  call  D  can  be  seen  to  contain  the 
absolute  values  of  the  eigenvalues  for  signals  propagating  to  the  solution  point  from 
either  left  or  right.  Indeed,  the  simplified  approximation  to  equation  (2) 

D6qj  =  A*q  ’  ,  -  (A*  -/»-)?  -  -O’.",  (13) 

leads  to  an  operationally  explicit  implicit  relaxation  procedure23  that  is  uncondi¬ 
tionally  stable  either  as  a  computed  interior  patch  boundary  point  or  general  interior 
point  relation.  Here  n,  n  + 1  means  data  from  either  time  level.  If  the  interior  point 
implicit  solution  procedure  is  two  level,  then  the  term  of  equation  (12)  at  the  inte¬ 
rior  point  j  —  1  or  j  + 1  will  be  linearized  (assumed  at  the  n  + 1  level)  as  in  equations 
(11a)  or  (lib). 

For  left  or  right  boundary  points,  the  frozen  (i.e.  not  computed  on  the  patch) 
data  at  j  —  1  or  j  +  1  in  equation  (13)  is  gotten  from  adjacent  patch  mesh  data. 
If  the  mesh  system  is  composite  and  mesh  lines  cross  the  boundary  to  the  solution 


point,  then  the  frozen  boundary  data  is  the  solution  point  data  of  the  adjoining 
mesl\. ,  For  the  case  of  nodes  on  lines  ending  at  the  patch  boundary,  which  relates 
equally  to  composite  mesh  with  lines  of  either  patch  ending  at  the  boundary  or  to 
overset  meshes,  the  frozen  data  is  got  by  interpolation  of  adjacent  mesh  patch  data 
to  the  boundary  point  location. 

Linear  interpolation  to  an  included  point  on  a  coordinate  line  or  within  a 
polygonal  cell  involves  only  positive  weights  on  the  interpolant  data.  In  most  of 
the  numerical  experiments  made  to  date  with  overset  grids  we  have  employed  a 
bilinear  interpolation31  based  on  the  four  comer  points  of  the  overlapping  mesh  cell 
enclosing  the  frozen  boundary  point.  However,  with  less  data  processing  a  linear 
interpolation  involving  the  three  corner  points  on  the  including  triangle  (Figure 
25a)  generalizes  to  the  use  of  the  four  corner  points  of  the  enclosing  tetrahedron  in 
three  dimensions. 

In  the  composite  mesh  case,  Figure  25b,  interpolation  is  naturally  along  an 
interior  coordinate  line  paralleling  the  zonal  boundary.  Such  interpolation  is 
one-dimensional  for  a  two-dimensional  problem  and  two-dimensional  for  a  three- 
dimensional  problem.  The  generic  composite  grid  problem  just  described  and  which 
we  have  tested  among  the  numerical  experiments  to  be  reported  in  the  next  section 
also  serves  as  a  gcdankcn  model  problem  for  the  overset  mesh  case.  Possible  solu¬ 
tions  that  come  to  mind  by  analogy  are  shown  in  Figure  25c  and  Figure  25d.  In 
both  figures  the  interpolation  is  along  two  point  coordinate  line  segments  (or  three 
point  triangular  surfaces  in  3-D)  of  the  adjacent  mesh  and  thus  is  a  direct  analog 
with  the  attendant  data  requirements  of  the  composite  mesh  case. 

We  use  stable  and  consistent  first  order  differencing  and  interpolation  at  the 
boundaries  regardless  of  the  order  of  accuracy  of  the  difference  approximation  in  the 
patch  interiors.  Since  the  divided  differences  of  the  computed  boundary  point  ap¬ 
proximation  are  of  the  same  accuracy  as  the  interpolation,  the  approaches  sketched 
in  Figures  25c  and  25d  may  be  regarded  as  letting  the  difference  operator  perform 
the  interpolation  (to  uniformly  spaced  data  Figure  25a)  in  the  direction  away  from 
the  computed  boundary.  Thus,  to  the  extent  that  the  solution  is  locally  well  rep¬ 
resented  by  a  linear  function  the  approaches  sketched  in  Figures  25a  and  25c  and 
25d  are  equivalent.  The  treatment  shown  in  Figure  25a,  however,  requires  the  same 
dimensionality  of  the  (triangular)  interpolation  procedure  as  that  of  the  problem 
and  one  higher  than  for  the  (linear)  treatment  of  Figures  25c  and  25d. 

As  a  final  theoretical  point  regarding  data  exchange  at  patch  mesh  boundaries, 
we  note  here  that  equations  (11)  imply30 

AF±  =  A±AF  (14) 

Hence  we  may  equivalently  write  the  right  hand  side  difference  operators  of  equation 
(13)  directly  in  terms  of  flux  components.  An  advantage  of  this  approach  is  that 
the  flux  components  normal  to  the  shock  discontinuity  are  continuous  across  the 


discontinuity.  Thus  interpolating  flux  data  from  one  grid  to  serve  as  needed  bound¬ 
ary  data  of  another  can  be  a  smoother  more  accurate  procedure  than  interpolating 
conservative  variable  data. 

We  close  this  section  by  noting  that,  consistent  with  the  2-D  interior  point 
schemes,17,33  we  difference  along  the  computed  boundary  coordinate  line  with  op¬ 
erators  written  for  the  boundary  aligned  coordinate.  In  the  diagonally  dominant 
approximate  factorizations  that  we  employ  in  multidimensions,  the  convective  mar 
trices  with  absolute  values  of  the  eigenvalues  for  both  coordinate  directions  are 
retained  in  the  diagonal  block. 

Two  Dimensional  Numerical  Experiments 

To  achieve  both  accuracy  and  robust  stability  in  this  difficult  problem  area, 
particularly  with  overset  meshes,  theory  can  provide  insight  into  what  to  attempt, 
but  the  acid  test  of  numerical  methods  is  performance  in  relevant  numerical  exper¬ 
iments. 

We  present  here  some  sample  results  with  major  findings  from  a  substantial 
number  of  experiments28  designed  to  test  various  aspects  of  the  accuracy  and  sta¬ 
bility  question  for  the  conservative  system  of  equations  for  gasdynamics  solved  on 
patched  meshes. 

The  numerical  experiments  to  be  discussed  here  involve  solution  of  an  inviscid 
flow  in  a  Mach  5  inlet  with  10°  compression  ramp  that  we  have  employed  in  previous 
experiments  with  first  and  second  order  upwind  methods  on  uniform  meshes17.  The 
problem  involves  two  of  the  generic  kinds  of  flow  structure,  shock  and  expansion 
fan,  which  are  not  possible  to  resolve  both  efficiently  and  to  the  extent  desired  on 
uniform  mesh.  As  the  result  of  interaction  of  the  expansion  fan  with  the  compression 
corner  and  reflected  shocks,  they  curve  in  non  simple  regions  for  which  the  exact 
solution  is  not  known  analytically.  The  coarse  base  level  grid  of  the  experiments 
has  26  x  26  points. 

Composite  Grid  with  Boundary  Overlap 

As  a  simple  test  of  employing  frozen  interpolated  boundary  data,  we  show  in 
Figure  26a  and  26b  the  grid  and  density  contours  for  a  patched  mesh  with  two 
full  cell  overlap  and  refinement  with  twice  as  many  mesh  points  in  the  streamwise 
direction  in  the  lower  patch.  Thus  every  other  mesh  point  at  the  upper  interior 
boundary  of  the  lower  patch  is  interpolated  between  computed  conservative  variable 
data  of  the  upper  patch  along  the  common  streamwise  coordinate  line.  This  test 
case  solving  sequentially  on  the  patches  with  the  first  order  method  is  numerically 
stable  with  local  time  steps  based  on  CFL  100.  There  is  no  oscillation  in  the  solution 
in  the  vicinity  of  the  patch  boundary  (shown)  and  refinement  in  the  lower  patch 
has  served  to  sharpen  the  solution  in  that  region,  though  high  gradient  regions  of 
the  solution  are  very  smeared  on  the  coarse,  nonaligned  meshes. 


Solutions  on  Uniformly  Refined  Mesh 

As  a  standard  for  comparison  with  results  from  overset  patched  refinement,  we 
show  in  Figures  27a  and  27b  pressure  contours  for  first  and  second  order  upwind 
methods  on  101  x  101  point  uniform  grids,  i.e.  4  times  refined  in  each  coordinate 
direction. 

Adaptive  Refinement  in  Overset  Grids 

Based  on  a  uniform  coarse  mesh  solution  similar  to  Figure  26b,  the  coarse  mesh 
was  overlaid  with  two  refined  mesh  patches  (Figure  28a)  aligned  with  the  compres¬ 
sion  comer  and  reflected  shocks.  Note  in  the  reflection  region,  the  two  overset 
patches  have  been  constructed  to  share  the  same  coordinate  lines  for  superior  grid 
communication.  The  coarse  grid  is  segmented  (broken)  under  the  overset  patches 
and  the  refinement  is  segmented  to  terminate  at  the  reflecting  boundary  (symmetry 
plane).  In  Figures  28b  and  28c  respectively,  we  show  pressure  contours  for  the  first 
and  second  order  upwind  methods  on  the  overset  grid. 

Discussion 

While  the  adaptive  refinement  about  doubled  the  coarse  mesh,  the  shock  struc¬ 
tures  treated  are  better  resolved  than  with  the  uniformly  refined  mesh  in  16  times 
the  points.  Thus  the  results  demonstrate  an  order  of  magnitude  improvement  in 
data  efficiency  to  be  gained  by  overset  refinement.  We  sketch  in  Figure  28d  a 
strategy  of  refinement  for  the  as  yet  untreated  expansions  of  the  problem. 

Graphics  for  Patched  Grids 

Graphics  is  an  important  tool  to  develop  and  debug  numerical  codes  and  an¬ 
alyze  numerical  results.  In  some  problems,  contour  plots  of  certain  flow  quantities 
such  as  pressure  and  Mach  number  are  sufficient  to  look  at;  in  other  problems,  one 
may  need  to  look  at  velocity  vectors  to  understand  the  solution  better.  In  the  case 
of  multiple  grids,  complication  arises  due  to  the  fact  that  solutions  in  more  than 
one  grid  are  available  in  some  regions  of  the  flow  domain  and  the  global  solution 
needs  to  be  composed  from  solutions  on  individual  grids. 

Graphics  for  Single  Grid  Solution 

Before  we  deal  with  multiple  grid  graphics,  we  outline  how  we  analyze  the 
single  grid  solutions.  Since  multiple  grid  solutions  are  made  up  of  single  grid  solu¬ 
tions,  insight  into  the  single  grid  graphics  will  help  to  understand  the  multiple  grid 
graphics  plotting  strategy.  In  the  process  of  solving  problems  there  are  three  roles 
for  graphics.  First,  the  grid  and  the  starting  solution  can  be  checked.  Second,  as 
the  solution  evolves,  graphics  can  be  used  as  a  tool  to  debug  and  understand  the 
time  evolution  of  the  solution  process.  Third,  the  converged  or  the  final  solution  is 
plotted  to  study  the  physics  and  compare  the  present  solution  with  other  solutions. 
Most  of  the  above  tasks  can  be  grouped  into  three  categories  of  plotting:  1)  Grid 
plotting;  2)  Contour  plotting  of  various  scalar  flow  quantities  and  velocity  vector 
plotting;  3)  Comparison  plotting.  Here  we  deal  with  the  first  two. 


Grid  Plotting 

Single  grid  plotting  is  done  by  essentially  drawing  straight  line  segments.  We 
use  the  commercial  macro  plotting  package  DISSPLA*  for  all  plotting  purposes. 
Our  plot  program  links  the  macro  DISSPLA  program.  Only  higher  level  commands 
need  to  be  defined  on  the  plotting  program  and  all  the  lower  level  commands  are 
defined  in  the  DISSPLA  program  and  are  transparent  to  the  users.  Given  the 
ordered  set  of  grid  points,  the  grid  plotting  routine  calls  the  line  segment  drawing 
command  for  every  line  segment  and  plots  the  grid.  Since  only  discrete  grid  values 
are  used  in  all  the  finite  difference  calculations,  the  representation  of  the  grid  by 
linear  elements  is  most  appropriate  and  no  interpolation  or  smoothing  is  necessary. 

The  same  philosophy  is  adhered  to  in  plotting  contour  curves  of  flow  quantities. 
Independent  of  the  order  of  the  scheme  and  the  solution  accuracy,  the  solution  is 
available  only  at  discrete  node  points  and  any  attempt  to  represent  the  solution  in 
a  smoother  sense  can  only  corrupt  the  solution.  With  this  in  mind,  we  use  a  fast, 
simple  and  robust  contour  plotting  routine. 

Contour  Plotting 

Contour  plots  for  the  two  dimensional  problems  are  very  useful  to  show  how 
accurately  shock  and  other  flow  structures  have  been  captured.  Any  scalar  flow 
quantity,  such  as  pressure,  density,  Mach  number,  temperature,  stream  function  or 
vorticity  can  be  plotted  for  various  constant  contour  levels.  The  general  method 
for  plotting  contour  levels  are  given  in  the  following  section.  First,  the  desired  flow 
quantities  are  calculated  from  the  dependent  variable  solution  vector  at  all  grid 
points  and  the  contour  subroutine  is  called  with  the  set  of  contour  values.  The 
contour  subroutine  computes  and  plots  the  various  contour  curves,  usually  in  the 
physical  domain. 

Since  any  finite  difference/finite  volume  formulation  solves  the  flow  field  in  an 
ordered  set  of  grid  points/cells,  the  construction  and  the  execution  of  the  contour 
program  is  structured  on  the  basic  grid  cell.  We  draw  all  contour  curves  cell  by  cell 
as  we  sweep  through  the  complete  grid.  At  the  corners  of  a  cell,  flow  function  F  has 
values  F 1,  F 2,  F3,  FA.  It  is  desired  that  the  contour  curve  for  function  level  FC 
need  to  be  plotted.  Then  along  each  side  of  the  base  grid  cell,  where  cross  over  points 
of  the  contour  curve  FC  occur,  they  can  be  found  through  linear  interpolation. 
By  connecting  sequentially  the  cross-over  points  encountered  among  the  sides,  one 
obtains  the  part  of  the  contour  curve  or  curves  in  a  given  cell  corresponding  to 
the  contour  value  FC.  By  repeating  this  process  for  all  cells,  the  complete  set  of 
contour  curves  for  the  whole  domain  can  be  found. 

Since  the  interpolation  used  to  find  the  contour  curve  in  the  base  cell  is  only 
linear,  the  accuracy  of  the  contour  curve  inside  the  cells  is  also  linear.  This  does 


*  DISSPLA  is  a  trademark  of  ISSCO. 


not  mean  that  the  solution  represented  by  the  contour  curves  is  first  order.  The 
discrete  solution  at  the  nodes  is  at  as  high  an  order  as  the  solution  technique.  If 
one  used  higher  order  interpolation  to  represent  the  discrete  data,  then  the  contour 
curves  represent  not  just  the  solution  but  non-physical/extra  smoothing.  Higher  or¬ 
der  interpolations  to  represent  discrete  data  can  also  result  in  an  oscillatory  solution 
representation  when  such  is  not  the  case. 

Velocity  Vector  Plotting 

Apart  from  the  contours,  at  times  it  is  also  desirable  or  necessary  to  look  at 
velocity  vectors.  This  may  be  to  study  the  location  of  separation  and  reattachment 
points  and  to  see  the  size  of  the  separation  zone.  The  development  of  the  bound¬ 
ary  layer  and  shear  layer  can  also  best  be  shown  by  velocity  vector  plots.  The 
present  graphics  code  provides  the  option  of  plotting  velocity  vectors.  To  plot  the 
velocity  vectors,  at  every  grid  point,  a  line  vector  with  or  without  an  arrow  head 
proportional  in  length  to  the  absolute  value  of  the  velocity  at  the  given  grid  point 
is  drawn.  Instead  of  the  velocity  vectors,  an  option  is  provided  to  plot  just  the 
velocity  direction.  In  problems  where  the  magnitude  of  the  velocities  may  change 
drastically,  it  has  been  found  convenient  to  plot  the  velocity  direction. 

Graphics  for  the  Multiple  Grids 

The  multiple  grid  graphics  code  is  based  on  the  single  grid  graphics  code.  In  the 
multiple  patched  grid  solution  procedure  grids  are  constructed  to  overlap.  When 
more  than  one  grid  is  used  to  solve  the  flow  problem,  the  major  question  that  arises 
is  what  to  do  in  regions  where  the  grids  overlap.  The  solution  is  obtained  in  all 
of  the  grids  and  so  there  are  regions  in  the  flow  field  where  multiple  solutions  are 
available.  The  accuracy  of  the  solution  on  each  grid  is  influenced  by  grid  fineness  and 
grid  topology,  among  overlapping  grids  the  solution  accuracy  can  be  quite  different. 
Since  the  solutions  in  overlapping  regions  do  not  have  the  same  function  values  or 
accuracy,  any  attempt  to  represent  them  there  will  exhibit  some  non-smoothness. 

Though  it  is  true  that  the  solution  in  the  overlapping  region  is  multivalued,  if 
the  solution  procedure  assures  smoothness  and  continuity  to  the  interface  boundary, 
then  any  one  of  the  solutions  in  the  overlapping  region  could  be  used  to  represent 
the  global  solution.  In  the  solution  procedure,  we  order  the  grids  with  assigned 
indices  in  some  sequence.  The  graphics  plotting  is  done  in  the  reverse  sequence.  In 
general,  the  coarsest  grid  is  the  first  grid,  and  any  finer  grid  interior  to  it  will  have 
a  higher  index  value.  To  plot  the  solution  in  the  global  domain,  the  grid  with  a 
highest  number  (finest  grid)  and  most  accurate  solution  will  be  plotted  first.  Next 
the  solution  in  the  lower  grids  will  be  plotted  sequentially.  In  regions  of  overlap, 
only  the  finer  grid  with  higher  index  values  is  used.  It  is  of  interest  to  note  that 
the  base  grid  is  subdivided  by  the  higher  level  grids  and  solution  in  the  finer  grid 
regions  (except  for  minimal  necessary  overlap)  are  available  only  from  the  finer  grid 
solutions  and  not  from  the  coarse  grid. 

In  Figure  28a  we  have  shown  an  example  of  overlapping  grids  for  a  supersonic 


inlet  problem.  Figure  28b  shows  the  pressure  contour  in  the  global  domain.  To  i 

obtain  this  plot,  the  pressure  contours  in  the  refined  and  shock  aligned  grid  were 

plotted  first  (with  the  single  grid  contour  plot).  Before  the  contours  in  the  base  grid 

were  plotted,  the  previously  plotted  refined  grid  region  was  “blanked  out”  so  that 

no  contour  lines  of  the  coarse  grid  solution  could  subsequently  be  plotted  inside 

the  region.  Blanking  a  curve  bounded  region  is  accomplished  without  special  effort 

using  a  feature  of  the  DISSPLA  graphics  program.  Finally,  the  base  grid  solution 

is  plotted  only  in  the  non-blanked  regions.  For  all  partitioned  grids,  like  the  base 

grid  in  this  case,  we  have  a  special  subroutine  that  plots  the  contour  curves  on  the 

integrable  cells  of  the  base  grid.  The  part  of  a  contour  line  outside  of  a  blanked 

region  for  a  cell  partially  overlapping  such  a  region  is  plotted  up  to  the  blanking 

boundary. 

The  above  choice  of  solution  representation  does  not  guarantee  smoothness 
and  continuity  of  contours  along  the  boundaries  of  overlapping  grids.  The  available 
discrete  solution  itself  is  nonsmooth  due  to  the  spatial  discretization  error  being 
different  in  each  grid.  Without  adding  additional  smoothing,  it  is  not  possible 
to  represent  the  multiple  grid  solution  with  smoothness  and  continuity.  We  have 
chosen  not  to  add  any  smoothing  but  to  represent  the  best  available  solution.  The 
technique  has  a  virtue  of  exhibiting  the  extent  of  truncation  errors  between  grids 
and  making  evident  places  where  more  refinement  may  be  needed.  In  our  numerical 
studies,  the  above  choice  of  contour  plotting  in  multiple  grids  seems  to  represent 
the  solution  very  well. 

Implicit  Upwind  TVD  Schemes  with  DDADI  Approximate  Factorization 

Methods  of  higher  than  first  order  accuracy  ii:  the  spatial  difference  approx¬ 
imations  are  not  monotonicity  preserving,  i.e.  lead  to  oscillation  and  nonlinear 
instability,  particularly  in  the  vicinity  of  shock  and  contact  discontinuities.  This 
problem  can  be  overcome  by  a  local  reduction  to  first  order  based  on  a  change  in 
eigenvalue  sign15,17  or  with  greater  generality  (and  effectiveness  in  multidimensions) 
by  the  use  of  flux  limiters.  Locally  one  dimensionally  monotonicity  preserving  or 
total  variation  diminishing  (TVD)  limiters  have  been  given  by  van  Leer32,33  for 
the  scalar  equation  with  Fromm’s  method.  The  first  of  these,  now  operating  on 
conservative  split  partial  flux  differences,  has  been  extended  to  systems  of  equa¬ 
tions  by  Yang34,35.  A  second  limiter  has  been  given  by  Harten3G  for  an  (unlimited) 

Lax-Wendroff  central  difference  method.  Harten’s  limiter  switches  to  upwind  dif¬ 
ferencing  in  expansions.  Similar  to  the  limiter  of  reference  33  for  Fromm’s  scheme, 
but  written  for  the  scalar  characteristics  representation  embedded  in  the  splitting 
of  systems  of  equations,  are  limiters  given  by  Chakravarthy  and  Osher37.  They 
have  given,  in  a  simple  one  parameter  format,  the  general  class  of  spatially  second 
order  biased  upwind  schemes  that  are  linear  combinations  of  central  and  upwind 
differencing.  The  above  three  limiters  were  extended  in  partial  flux  difference  split 
pieces  to  systems  of  equations  and  adapted  to  implicit  methods  in  reference  38. 
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Unconditionally  stable  linearized  block  implicit  schemes  in  delta  form  can  be 
constructed,  e.g.  Warming  and  Beam30,  with  first  order  upwind  differencing  on  the 
left  hand  side  (LHS)  of  the  difference  equations  and  second  order  biased  upwind 
(flux  limited)  differencing  on  the  right  hand  side  (RHS).  Approximate  factorizations 
of  the  LHS  block  matrix  lead  to  labor  and  storage  efficient,  locally  one  dimensional 
solution  procedures  that  may  be  block  tridiagonal,  block  bidiagonal,  or  approxi¬ 
mately  scalar17,40,41. 

Implicit  schemes  with  LU  decompositions  for  hyperbolic  equations  were  con¬ 
structed  by  Jameson  and  Turkel42.  One  of  their  important  findings  is  that  the 
resulting  three-dimensional  alternating  direction  algorithm  is  unconditionally  sta¬ 
ble  in  the  linear  case  and  also  yields  a  steady-state  solution  which  is  independent 
of  time  step  so  long  as  each  factor  of  the  L  and  U  decompositions  is  diagonally 
dominant. 

For  upwind  schemes,  Lombard,  et.  al.17  have  given  a  simple  diagonally  dom¬ 
inant  approximate  factorization  (DDADI)  that  retains  the  contributions  (absolute 
values  of  the  eigenvalues)  for  all  coordinate  directions  on  the  diagonal  block  in  each 
factor.  The  approach  leads  to  ADI  block  tridiagonal  and  block  bidiagonal  schemes 
that  have  been  found  in  numerical  experiments  to  be  more  robust  than  approxi¬ 
mate  factorizations7,43  in  which  eigenvalue  contributions  are  separated  among  the 
factors.  The  DDADI  approximate  factorization  has  also  led  to  a  class  of  signifi¬ 
cantly  more  rapidly  convergent  single  data  level  algorithms23  that  use  symmetric 
Gauss-Seidel  relaxation. 

In  the  present  study,  following  the  DDADI  approach,  we  describe  new  upwind 
second-order  TVD  implicit  algorithms  for  the  Euler  and  Navier-Stokes  equations. 

Some  preliminary  computational  results  are  presented  for  2-D  and  axisymmet- 
ric  flows. 

Numerical  Algorithms 

We  consider  the  2-D  Euler  equations  of  gasdynamics  in  general  curvilinear 
coordinate  system 

drQ  +  d^F  d^G  =  0  (15) 

where  Q  =  Q/J  and  F  =  {ZxF+tyG)/J,  G  =  (r]xF-\ri]yG)/J,  and  J  =  CxVv-^x, 
the  metric  Jacobian.  Where  Q  =  (p,  pu,pv,e)T  is  the  conservative  variables,  F  = 
(pu,pu2  +  p,puv,u[e  +  p))T  and  G  =  [pv,  pvu,  pv2  +  p,  v[e  +  p))T  are  the  flux 
vectors  in  Cartesian  coordinates.  Here  p  is  the  gas  density,  u,  v  are  the  gas  velocity 
components,  p  the  pressure  and  the  total  energy  e  =  p/( 7  —  1)  -I-  \p{u2  +  v2). 

We  employ  splittings  of  the  adjacent  two  point  conservative  flux  difference  in 
the  Lombard- Yang  generalized  Roe  forms44,45 

A  F±  =  A±AAQ  =  A±AQ  (16a) 

=  A±AF  =  A±AAQ  (166) 
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with  the  projection  operators  A±  =  ( M  T  I*  T  1  M  *)  providing  partitioning 
according  to  the  signs  of  the  eigenvalues  by  =  |(/±  sgn(A)).  Here  the  matrices 

M  1  and  T  1  give  the  discrete  difference  transformations  between  conservative  and 
primitive,  and  primitive  and  characteristic  variable  representations  respectively. 

With  two  splittings  CSCM17  and  CFDS35  we  can  construct  implicit  upwind 
schemes,  respectively,  linearizing  the  flux  differences  in  the  forms  (16a)  and  (16b). 

The  first-order  upwind  implicit  scheme  using  backward  Euler  time  differencing 
and  linearization  based  on  the  form  16b,  for  example,  can  be  given  as 


\I  +  +  kA7+t  **i*i*A  + 

=  RHSj  +  RHS,  =  RHS  (17) 


where 


and 


RHS{  =  -hA^A^F  -  hA-Hi  kAi+i,kF 
RHS,  =  -hB+k_^AJk_kG  -  hBjk+kAf  k+iG 


Here  h  is  the  time  step,  and 


«<??> = -  Qj* 

A j-k,kF  =  Fy.fc  -  A =  Gj.k+1  -  Gj,k 

Block  Tridiagonal  Approximate  Factorization  for  the  Implicit  LHS 

We  will  couch  the  remainder  of  the  discussion  of  some  efficient  solution  proce¬ 
dures  in  a  simplified  notation  based  on  the  form  16a  with  obvious  correspondence 
to  the  form  16b.  Then  equations  (17)  can  be  expressed  in  the  matrix  form 

\-A+  ,  -B+  ;  I+\A\  +  \B\  ;  B~  ,  A~  \6Q  =  RHS  (18) 


The  LHS  of  equation  (18)  may  be  approximately  factored  in  the  DDADI  form 

[  ~A+  ,  D  ,  A-  ] D-l[-B+  ,  D  ,  B~  \6Q  (19) 

with  D  =  I  +  \  A\  +  \B\.  Equation  (19)  gives  the  block  tridiagonai  solution  sequence 

[  -A+  ,  D  ,  A~  1 6Q*  =  RHS  ,  x 

(  -B+  ,  D  ,  B~  ]<5 Qn  =  D6Q* 

Qn+ 1  =  Qn  +  6Qn 
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Second  Order  TVD  Schemes 

The  RHS  of  the  above  schemes  based  on  equations  (17)  is  only  first-order  in 
time  and  space,  and  in  general  is  not  sufficient  for  flow  resolution.  In  this  study 
we  further  investigate  the  possibility  of  replacing  the  RHS  by  robust  second-order 
methods  and  retaining  the  unconditionally  stable  LHS  to  achieve  an  efficient  algo¬ 
rithm  for  steady-state  calculations. 

In  the  following  we  give  two  second-order  flux  limited  difference  approximations 
written  for  the  £  direction.  As  is  common  practice,  the  schemes  are  written  as  a 
first  order  method  with  added  flux  limited  correction  terms  that  in  the  absence  of 
limiting  render  the  method  second-order  accurate. 

Scheme  I,  Adapted  Lax-Wendroff/Harten 

First,  a  modified  flux  vector  is  defined  for  a  hyperbolic  system  of  conservation 
laws  based  on  the  characteristic  flux  difference  splitting  concept.  Then  the  limited 
second  order  method  couched  as  a  first  order  upwind  method  is 

RHS{  =  ~hA;,i  kAM,kFu  -  hA-+i  tAj+i,kF“  (21) 

where  FM  is  the  modified  flux  vectors.  The  value  of  FM  at  nodal  point  j,  k  is  given 
by 

Ffi  =  Fjj.  +  (22) 

where  F  is  the  original  flux  vector  and  E  is  an  additional  vector  of  second  order 
correction  terms  remaining  to  be  defined. 

The  column  vector  E  at  nodal  point  j,k  is  Ejt^  =  (c1>  fc ,  C2j  k , ...,  *Aj<k)T  with 
its  /  components  given  by 

'</..  =  s'Jt4,.min(l?',+4..M'''j-i.J)-  when  -  0> 

=  0,  when  0,  (23) 

where  <?/  .  (/  =  1,2, ...,  4)  are  components  of  the  following  column  vector 

JT^  i» 

-%’+$,*  =  s^-4i+$,fcAy+$)*ir/2  ~  kAJ+.£)fcAJ+£ifcF/2  (24) 

and 

‘w. = sgn(iW.>  <25> 

The  sgn  A  in  equation  (24)  is  given  by 

sgn  o  =  MT{diag{sgn  (26) 
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In  equation  (21)  with  (22)  the  first  order  terms 

A±AF  =  A  F±  (27) 

can  be  represented  alternatively  using  the  flux  vector  splittings  of  Steger  and 
Warming43  or  van  Leer46.  The  latter  introduces  numerical  dissipation  that  ren¬ 
der  solutions  in  the  vicinity  of  shocks  and  contacts  monotone. 

Scheme  I  can  be  implemented  as  an  implicit  relaxation  scheme  as  in  equation 
(17)  or  as  a  second  order  time  accurate  time  split  explicit  method  given  by 

=  LtL„L„L(Qlk  (28) 

For  the  operator  in  the  ^-direction,  we  have  for  example  the  first  intermediate 
result 

<3;,*  =  G?,*  +  RHS{  (29) 

with  RHS(  given  by  equation  (21).  Results  for  Scheme  I  in  an  implicit  relaxation 
method  applied  in  transonic  flow  have  been  given  in  reference  38.  Here  we  give  re¬ 
sults  for  the  explicit  scheme  with  van  Leer46  flux  splitting  for  the  first  order  upwind 
terms47  applied48  to  2-D  unsteady  strong  shock  diffraction  around  a  hypersonic  or¬ 
bital  transfer  vehicle  (AOTV)  configuration.  The  unsteady  problem  has  potential 
relevance  both  in  maneuver  and  in  simulating  the  effect  of  nonuniformity  in  the 
upper  atmospheric  flight  regime. 

Figure  29  shows  the  body  shape  in  a  computational  grid  generated  with  the 
interactive  algebraic  procedure.  In  Figures  30a  and  30b  respectively  we  show  pres¬ 
sure  and  Mach  contour  plots  at  an  instant  in  time  for  a  case  rim  at  Mach  30,  7  = 
1.1  and  15°  angle  of  attack.  Other  work49  has  shown  important  dependencies  such 
as  shock  standoff  distance  and  temperature  on  the  value  of  7.  In  future  work  we 
will  introduce  aspects  of  the  present  accurate  and  robust  numerical  methodology 
into  upwind  codes  incorporating  variable  7  based  in  equilibrium  and  nonequilibrium 
reacting  gas  procedures. 

Scheme  II,  Adapted  and  Extended  Fromm/van  Leer,  Chakravarthy-Osher 
The  scheme  is 

(RHS)s  =  -hAFt^  -  hAF-+i 

+  -  A F*+i)  +  (A Fj+i  -  A?;_4)  (30) 

where  the  A F  and  A F  are  respectively  backward  and  forward  limited  flux  differ¬ 
ences.  Following  reference  37  for  example, 

Afj+J  =  (M  f)ft,mmmod  (A/++4,3A/t_4) 
and 

A F*_t  =  (MT)S_ iminmod  (A/+_4,3A/t_4). 


(31) 


Here  the  characteristic  partial  flux  differences  are 


A /*  =  I±  T  1  M  1  A F 


(32) 


and 

minmod(:r,  y)  =  sign(x)  max(0,min(|x|,  y  sign(x))]. 

As  an  alternative  to  limiting  the  characteristic  partial  flux  differences  we  also  con¬ 
sider  here  a  new  scheme  limiting  in  the  conservative  split  partial  flux  differences,  in 
conformity  with  scheme  I.  Then  we  replace  equation  (31)  with,  for  example 

AF^"+ j  =  minmod(AF^^,3Ai?y^i)  (33) 

Scheme  I  can  be  disected  to  show  that  in  regions  of  compression  the  scheme 
yields  a  second  order  central  difference  approximation  and  in  regions  of  expansion, 
second  order  upwind.  Chakravarthy  and  Osher37  gave  a  one  parameter  family  of 
flux  limited  second  order  schemes  that  can  be  seen  to  be  linear  combinations  of  the 
central  and  upwind  methods.  The  Fromm  scheme  given  here  for  illustration  is  the 
case  of  ^  central  and  ^  upwind  methods.  Another  case  with  |  central  differencing 
and  |  upwind,  yields  a  third  order  accurate  method. 

Computational  experiments  with  the  one  parameter  family  of  second  order 
schemes  for  the  advective  terms  of  the  compressible  Navier-Stokes  equations  have 
been  performed  in  the  context  of  a  much  tested  45°  -  15°  transonic  nozzle  problem, 
Figure  31.  Figure  32a  shows  Mach  contours  in  the  vicinity  of  the  throat  for  the 
unlimited  second  order  upwind  method.  Limiting  in  both  the  characteristic  and 
conservative  partial  flux  difference  pieces  has  been  carried  out.  Here  we  show  in 
Figures  32b  and  32c  Mach  contours  for  the  Fromm  scheme  with  the  two  approaches 
to  limiting.  The  unlimited  upwind  scheme  captures  the  weak  shock  somewhat  more 
sharply  but  is  considerably  less  robust.  Limiting  in  the  conservative  flux  split  pieces 
is  simpler  to  implement  and,  taking  into  account  eigenvector  variation,  may  have 
better  nonlinear  stability  than  limiting  in  the  characteristic  representation;  but  the 
solution  also  appears  less  smooth  behind  the  weak  shock. 

Application  of  3-D  Single  Level  Algorithm 
to  Multiple  Nozzle  Flow  on  Segmented  Composite  Mesh 

A  clustered  multi  rocket  nozzle  shrouded  exhaust  system  provides  a  problem  of 
relevant  kind  and  geometric  complexity  typical  of  those  needing  advanced  compu¬ 
tational  techniques  to  analyze.  Here  the  symmetric  Gauss-Seidel  implicit  method 
of  planes  algorithm  described  earlier  is  adapted  to  a  segmented  composite  grid. 

The  grid  covers  a  “pie  slice”  sector  bounded  by  symmetry  planes  through  a 
circular  array  of  engines  around  a  single  engine  at  the  center.  The  symmetry  planes 
pass  through  the  axis  of  the  center  engine  and  either  through  or  between  the  nozzles 
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of  the  outer  engines.  Figure  33a  shows  a  3-D  perspective  view  of  the  surface  mesh 
of  the  computational  domain.  Figure  33b  shows  a  symmetry  plane  mesh  through 
the  center  and  an  outer  engine,  and  Figure  33c,  the  projection  of  a  cross  stream 
mesh  surface  in  the  shrouded  exhaust  region.  Figure  34,  which  relates  to  Figure 
33c,  shows  the  segmented  composite  mesh  coordinate  topology.  Where  the  index  L 
is  for  the  azimuthal  coordinate  0  of  the  outer  engine  associated  mesh  patch  and, 
similarly,  K  is  for  the  radial  like  coordinate  rjt  the  center  engine  associated  mesh 
patch  sector  is  naturally  mapped  onto  a  segment  of  the  outer  engine  grid  topology 
with  continuity  of  the  rj  coordinate  in  K 1  +  1  <  K  <  KM  AX  and  1  <  L  <  Ll. 

Boundary  Conditions 

The  segmented  composite  mesh  approach  requires  the  facility  to  impose  differ¬ 
ent  boundary  condition  procedures  on  limited  segments  of  coordinate  curves  piece- 
wise  fitted  to  boundaries.  Segments  of  viscous  wall  boundary  are  treated  with 
implicit  characteristic  based  procedures  as  described  in  reference  17.  Boundary 
segments  of  inviscid  symmetry  planes  can  be  treated  by  differencing  to  out  of  do¬ 
main  data  obtained  by  applying  symmetry  relations  operationally  and  explicitly 
but  with  unconditional  stability  based  on  DDADI  as  explained  earlier  in  this  report 
for  interior  patched  mesh  boundaries.  Alternatively,  symmetry  boundaries  can  be 
treated  completely  locally  like  inviscid  walls  by  an  implicit  matrix  procedure  involv¬ 
ing  characteristic  relations  to  the  boundaries,  augmented  by  local  invariance  of  the 
zero  normal  velocity  component  as  in  reference  17,  and  extrapolations  of  tempera¬ 
ture  and  tangential  velocity  components  through  constructed  stable  characteristic 
like  relations  as  described  in  reference  50. 

Singular  poles  such  as  at  K  =  1  and  KM  AX  of  the  grid  topology  are  buried  at 
half  cell  spacing  in  the  mesh.  Operationally  explicit  DDADI  differencing  is  carried 
out  from  first  interior  points  to  cross  pole  data  obtained  from  symmetry  conditions 
as  given  in  Figure  34  for  the  pole  at  K  =  1. 

Figure  35  shows  Mach  contours  on  a  cross  flow  coordinate  surface  intersecting 
the  backward  facing  step  outside  the  exit  of  the  outboard  nozzle.  This  region  is 
impacted  by  the  complicated  interaction  of  a  weak  shock  exiting  the  nozzle  and 
an  expansion  fan  whose  strength  at  the  nozzle  exit  corner  is  influenced  by  the 
size  of  the  local  backward  facing  step.  Figures  36a,  b  and  c  show  Mach  contours 
for  longitudinal  mesh  surfaces  associated  with  the  outer  nozzle  flow.  The  surfaces 
intersect  the  inter  nozzle  symmetry  plane,  respectively,  in  the  Ll  corner,  at  the 
minimum  step  between  outboard  nozzles,  and  at  the  maximum  step  in  the  corner 
where  the  shroud  intersects.  The  latter  figures  show  the  nozzle  weak  shock,  the  base 
corner  expansion  and  then  the  recompression  shock  against  the  symmetry  plane, 
and  finally  the  extensive  flow  separations  in  the  corners.  The  extent  of  departure 
from  axisymmetry  in  the  plume  region  caused  by  the  azimuthally  varying  base  step 
heights  can  be  appreciated.  Figure  36d  shows  Mach  contours  in  the  symmetry  plane 
through  the  inner  and  outer  nozzle  centers.  Here  appears  the  turbulent  viscous  wake 


of  the  inter-nozzle  exit  step  as  another  captured  feature  of  the  multi-nozzle  plume 
flow. 

Computational  Experiments  in  Patched  Mesh  Topology 
and  Solution  Accuracy 

For  high  Reynolds  number  viscous  flows,  particularly  for  massively  separated 
base  flows  computed  in  the  efficient  Baldwin-Lomax  thin  layer  approximation,  mesh 
topology  has  a  great  deal  to  do  with  accuracy  of  the  resulting  solution.  Occurring 
in  numerous  aerospace  applications,  separated  base  flow  is  of  significant  interest  to 
predict  for  it’s  impact  on  drag  and  moments  on  aircraft,  projectiles  and  tactical 
missiles.  It  is  of  interest  for  efficient  design  of  dump  combustors  and  to  avoid 
chemical  erosion  problems  in  liquid  propellant  injector  systems.  Base  flow  is  also  of 
interest  for  thermal  protection  requirements  on  multi-rocket  powered  vehicles  such 
as  the  Space  Shuttle  and  hypersonic  reentry  bodies  such  as  the  Galileo  Probe  and 
the  Aeroassisted  Orbital  Transfer  Vehicle  (AOTV). 

In  recent  years,  the  evolving  computational  capability  for  viscous  flows  has 
led  to  attempts  to  analyze  such  base  flows  numerically.  The  widely  perceived 
need  to  assess  the  accuracy  of  the  techniques  particularly  for  solving  the  compress¬ 
ible  Navier-Stokes  equations  for  turbulent  base  flows,  has  led  to  experiments  with 
broadly  participated  computational  efforts  featuring  axisymmetric  afterbody  flows 
with  a  centered  propulsive  jet.  The  time  dependent  Navier-Stokes  solutions  have 
been  computed  by  Deiwert,  Andrews  and  Nakahashi51  and  Sahu  and  Nietubicz52. 
A  novel  feature  of  that  work  is  the  use  of  flow  structure  adaptive  meshing. 

Of  particular  relevance  here  are  computations  of  Sahu  and  Nietubicz52  and 
Thomas,  Reklis,  Roloff  and  Conti53  for  the  MICOM  model  tactical  missile 
experiment54  at  freestream  Mach  number  1.4.  The  axisymmetric  flow  experiment 
featured  a  tangent  ogive  cylinder  body,  Figure  37,  terminating  in  a  square  cylindri¬ 
cal  base.  A  centered  nozzle  with  flush  exit  occupied  only  0.2  base  diameters  of  the 
model  which  was  9  diameters  long.  The  nozzle  exit  Mach  number  was  2.7  for  an 
underexpanded  jet  with  pressure  ratio  2.15  relative  to  freestream. 

Sahu  and  Nietubicz52  employed  the  Beam-Warming7  central  difference  algo¬ 
rithm  with  the  Baldwin-Lomax8  algebraic  turbulent  eddy  viscosity  model  and  the 
thin  layer  Navier-Stokes  approximations.  Since  viscous  terms  were  not  applied  in 
the  streamwise  direction,  the  base  wall  was  treated  inviscidly.  Thomas,  et  al.53,  who 
employed  the  Thomas  and  Lombard  version  of  the  central  difference  algorithm,  ap¬ 
plied  the  full  Navier-Stokes  equations  with  eddy  viscosity  computed  from  alternative 
two  equation  turbulence  models.  Both  groups  of  computations  were  made  without 
knowledge  of  the  experimental  results  except  for  data  given  along  a  line  1.6  diam¬ 
eters  upstream  of  the  base  to  serve  as  an  inflow  computational  boundary.  Unused 
data  was  also  given  on  a  lateral  boundary  1.6  base  radii  from  the  model  axis. 

Both  insight  and  computational  experiments  suggest  there  are  two  ways  that 
topology  influences  solution  accuracy.  The  first  is  through  the  efficiency  of  utiliza- 


tion  of  mesh  point  resources  to  achieve  balanced  resolution.  The  second  is  through 
alignment  with  and  natural  accurate  capture  of  flow  structures  on  the  mesh. 

Here  we  report  on  the  first  two  parts  of  a  continuing  capsule  study  of  the  impact 
of  computational  meshing  strategies  on  solution  accuracy  in  propulsion  base  flow. 
The  MICOM  model  problem  with  the  experimental  and  previous  computational 
data  base  was  chosen  as  the  vehicle  for  the  study.  The  first  part  of  the  study 
dealt  with  two  complementary  composite  mesh  topologies.  The  first  is  a  dual  wrap¬ 
around  composite  mesh  constructed  in  two  blocks.  The  second  is  a  composite  step 
mesh  constructed  in  three  blocks.  Based  on  the  techniques  described  earlier  in  this 
report,  the  blocks  of  smooth  mesh  are  generated  algebraically  in  subpatches. 

The  step  mesh,  shown  in  Figure  38b,  is  similar  to  those  employed  in  references 
52  and  53,  and  naturally  fits  and  is  simple  to  construct  for  (near)  right  angle  cor¬ 
ners  such  as  occur  at  the  cylinder  wall  and  nozzle  exit  junctures  with  the  base.  The 
wrap  around  mesh,  Figure  38a,  provides  a  continuity  of  the  thin  layer  viscous  ap¬ 
proximation  around  the  base  corners  and  is  thought  to  generally  make  most  sense 
of  locally  one  dimensional  algebraic  eddy  viscosity  modeling  (modified  Baldwin- 
Lomax)  which  we  employ  here.  Details  of  the  turbulence  model  and  the  application 
of  the  CSCM  method  are  given  in  reference  55. 


Propulsion  Base  Flow  Results 

For  the  two  base  flow  meshes,  we  present  results  in  Mach  number  contour  plots. 
The  results  for  the  two  meshes  in  comparison  with  each  other  and  with  the  results 
of  Sahu  and  Nietubicz52  and  Thomas,  et  al.63  are  quite  interesting  and  informative. 

In  Figures  39a  and  39b  we  show  Mach  contour  plots  for  the  wrap-around  and 
step  base  flow  meshes.  The  two  plots  exhibit  the  same  flow  structures  -  expansions 
at  the  nozzle  and  upper  base  corners,  followed  by  weak  turning  (recompression) 
shocks  in  the  outer  wake  flow  and  inner  jet  flows,  large  stagnant  separation  bubbles 
in  the  base  region,  Mach  disc  and  reflected  shock  in  the  jet  flow,  and  finally  wake 
expansion  downstream  of  the  recompression  neck  behind  the  turning  shocks. 

In  addition,  we  see  shear  layers  emanating  from  the  nozzle  and  upper  base 
corners  and  turning  with  the  recompressing  flow  to  be  carried  down  the  wake.  Also 
evident  in  the  figures  is  the  effects  of  the  structure  of  the  slow  recirculating  flows 
in  the  separated  base  region.  All  the  above  structure  is  in  very  good  qualitative 
agreement  with  that  shown  in  plots  in  reference  52  and  53. 

Considering  the  contour  plots  and  the  associated  meshes  among  our  results  and 
those  referenced,  two  observations  come  to  mind.  The  first  is  that  different  solutions 
obviously  express  different  structures  with  greater  or  lesser  accuracy.  The  second 
is  that  all  the  meshes  employed  here  and  in  the  references  cited  employ  locally 
one  dimensional  clustering  and  stretchings  that  emphasize  resolution  in  certain 
regions  or  directions  at  the  expense  of  others.  Where  grid  points  are  well  distributed 
relative  to  resolving  a  particular  flow  feature,  it  is  sharp  and  (presumably)  well 
located.  Where  not,  the  converse  holds.  Another  more  significant  mesh  dependent 
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feature  is  apparent  in  the  formation  region  of  the  shear  layers  from  the  nozzle  base 
corner.  In  the  wrap-around  mesh  case  of  Figure  38a,  the  boundary  layer  separates 
a  short  distance  around  the  corner  at  the  base  wall,  whereas,  in  the  step  mesh  case 
the  separation  appears  sharply  at  the  comer.  Associated  with  this  mesh  topology 
dependent  phenomena  are  observable  differences  in  the  angle  and  possibly  energy  of 
the  expansion  jet  down  the  near  base  shear  layer.  The  wrap-around  grid  has  much 
better  resolution  near  the  tail  of  the  separation  bubble,  and  also  generally  better 
resolution  over  shear  layers  on  either  side  of  the  bubble. 

What  one  might  term  simple  composite  mesh  schemes  of  the  type  employed  in 
Figures  38a  and  38b  have  the  obvious  flaw  that  they  are  comparatively  inflexible 
with  regard  to  continued  grid  point  clusterings  throughout  the  mesh  once  clustering 
has  been  initiated  of  necessity  somewhere,  as  for  example  in  boundary  layers.  The 
rigidity  leads  to  wasted  mesh  point  resources  in  some  inappropriately  clustered 
regions  at  the  expense  of  needed  resolution  in  others. 

The  techniques  of  multiple  patched  meshes  with  overlap  to  stably  exchange 
data  through  interpolation  at  interior  patch  boundaries  as  described  earlier  in  this 
report  offer  the  freedom  to  utilize  mixed  topologies  and  commit  mesh  resources 
as  appropriate  to  various  zones  of  the  problem.  The  second  part  of  the  study 
employed  the  multiple  patched  grid  techniques  to  effect  globally  a  hybrid  topology 
of  the  wrap-around  and  step  mesh  types. 

Figure  40a  shows  the  patched  grid  in  the  large  where  it  appears  as  a  step  mesh 
but  with  less  clustering  in  the  base  wake  region  vis  a  vis  Figure  38b.  The  wrap¬ 
around  feature  is  incorporated  in  patches  of  refined  overset  mesh  that  are  local  to 
the  upper  and  lower  exterior  comers  at  the  base.  Figure  40b  shows  a  local  detail  of 
the  overset  patch  boundary  at  the  lower  (nozzle  exit)  comer.  Also  shown  are  the 
overlapping  boundaries  of  other  grid  patches  that  intersect  in  the  vicinity. 

Figure  41  shows  a  Mach  contour  plot  for  the  solution  on  the  hybrid  grid.  By 
comparison  with  Figures  39a  and  39b,  all  flow  structure  features  of  the  solution  on 
the  multiple  patch  grid  appe  more  sharply  resolved  though  the  total  number  of 
points  in  each  global  grid  is  comparable.  Finally,  Figure  42  provides  a  comparison 
of  base  pressure  results  among  computations  and  the  experiment  of  Petrie  and 
Walker54.  Except  near  the  nozzle  exit  corner,  the  results  for  our  computation  on  the 
step  mesh  of  Figure  38b  is  quite  close  to  that  of  Sahu  and  Nietubicz52  and  to  reduce 
confusion  is  not  plotted.  All  the  computed  results  show  a  similar  wave  pattern  in  the 
base  associated  with  the  two  counter-rotating  circulation  cells.  That  this  structure 
is  not  apparent  in  the  experimental  results  raises  some  questions.  Further  study 
of  this  problem  as  to  mesh  topology,  adaptive  refinement  and  turbulence  modeling 
will  prove  useful  in  advancing  understanding  in  the  difficult  and  important  area  of 
base  flow  computation. 
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Dr.  N.  Nagaraj 

Dr.  J.  Y.  Yang 

Dr.  R.  C.-C.  Luh 

4.  Interactions 

The  research  described  in  this  report  has  to  this  point  been  partially  presented 
in  the  form  of  a  paper  on  algebraic  grid  generation  by  Vinokur  and  Lombard  (refer¬ 
ence  5),  a  SIAM  meeting  paper  (reference  19)  by  Oliger  and  Lombard  on  boundary 
procedures  for  bidiagonal  alternating  sweep  schemes  and  a  Computational  Fluid 
Dynamics  Seminar  at  NASA-Ames  Research  Center  by  Lombard  on  the  Univer¬ 
sal  Single  Level  Implicit  Algorithm.  The  latter  invoked  tremendous  interest  and 
discussion. 

The  research  has  spawned  four  meeting  papers  on  the  single  level  relaxation 
algorithm  -  references  23,  24  and  25  and 

Lombard,  C.K.,  Raiszadeh,  Farhad  and  Bardina,  Jorge:  “Efficient,  Vectoriz- 
able  Upwind  Implicit  Relaxation  Algorithms  for  Three-Dimensional  Gasdynamics,” 
SIAM  Spring  Meeting,  Pittsburg,  PA,  June  1985. 

Six  papers  have  been  written  on  interior  patch  boundary  treatment,  data  struc¬ 
ture  and  applications  of  patched  mesh  systems.  These  are  references  28,  29,  31  and 
55  of  the  present  report  and 

Lombard,  C.K.  and  Venkatapathy,  Ethiraj:  “Implicit  Boundary  Treatment  for 
Joined  and  Disjoint  Patched  Mesh  Systems,"  Workshop  on  High  Reynolds  Flows, 
Nobeyama,  Japan,  September  1985. 
and 

Lombard,  C.K.,  Bardina,  J.,  Venkatapathy,  E.,  Yang,  J.Y.,  Luh,  R.C.-C.,  Na¬ 
garaj,  N.  and  Raiszadeh,  F.:  “Accurate,  Efficient  and  Productive  Methodology  for 
Solving  Turbulent  Viscous  Flows  in  Complex  Geometry,”  Presented  at  the  10th  In¬ 
ternational  Conference  on  Numerical  Methods  in  Fluid  Dynamics,  Beijing,  China, 
June  1986  and  soon  to  be  published  in  a  bound  volume  Lecture  Notes  in  Physics. 

The  research  has  given  rise  to  two  papers  on  the  3-D  CSCM-S  relaxation  al¬ 
gorithm  with  application  to  multi  rocket  engine  exhaust  flow,  a  problem  of  interest 
to  and  partially  supported  by  AFRPL.  The  papers  are 


Bardina,  Jorge  and  Lombard,  C.K.:  “Three  Dimensional  CSCM  Method  for  the 
Compressible  Navier-Stokes  Equations  with  Application  to  a  Multi-Nozzle  Exhaust 
Flowfield,”  AIAA  85-1193,  June  1985. 
and 

Lombard,  C.K.,  Bardina,  Jorge  and  Raiszadeh,  Farhad:  “Shrouded  Multi- 
Rocket  Engine  Plume  Flow  by  an  Efficient  Upwind  Relaxation  Algorithm  for  the 
Compressible  Navier-Stokes  Equations,”  annual  meeting  of  the  JANNAF  Exhaust 
Plume  Technology  Subcommittee,  San  Antonio,  Texas,  May  1985. 

5.  New  Discoveries 

The  Universal  Single  Level  Algorithm  for  the  compressible  Euler  or  Navier- 
Stokes  equations  is  a  new  discovery  in  numerical  methods  that  promises  to  result 
in  substantial  efficiencies  in  data  storage,  programming,  machine  time  and  human 
productivity. 
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Figure  1.  Portion  of  an  algebraically  generated 
computational  mesh  for  a  flow  domain  containing  an 
external  and  an  internal  corner. 
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Figure  2.  Representative  2-D  mesh  generated  with 
simplified  interactive  algebraic  procedure. 


Figure  3.  Representative  3-D  mesh  generated  with  simplified 
interactive  algebraic  procedure. 


Figure  4.  Shubin's  diverging  nozzle  supersonic 
flow  solution  developed  in  one  forward  sweep  from 
supersonic  initial  data. 


Figure  5.  continued 


Figure  6.  Blottner's  converging-diverging  supercritical  nozzle  flow  solution 
developed  in  alternating  direction  sweeps  for  ten  global  iteration 
steps.  square  -  computed  results  ,  line  -  exact  solution 


Figure  6.  continued 
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Figure  7.  Subcritical  solution  to  Blottner's  c 
developed  with  alternating  sweeps  in 
square  -  computed  results  ,  line  - 


Figure  8.  Solution  to  subcritical  nozzle  problem  after 
15  global  iterations  with  forward  marching  sweeps  only. 


relx  =  0.010 


Figure  9.  Solution  to  subcritical  nozzle  problem  after 
6C  global  iterations  with  forward  marching  sweeps  only. 
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Figure  16.  Pressure  contours  from  the  first  order  viscous 
solution  for  the  inlet  after  20  global  sweeps  with  a  51  x  51 
stretched  mesh.  Note  the  leading  edge  shock. 
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Figure  17.  Wall  pressure  comparison.  Line  -  exact  inviscid 
solution,  square  -  first  order  inviscid  solution, 
dot  -  first  order  viscous  solution. 
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Figure  18.  Convergence  history  of  the  RMS  of 
the  residuals  for  the  inviscid  first  order  inlet 
problem  with  4  inner  iterations  per  global  sweep. 
Note  at  the  end  of  the  first  sweep  the  residual 
has  reduced  and  the  solution  has  converged  for  all 
Practical  purposes. 


Figure  19.  Perspective  view  of  the  wall  surface  mesh  for 
the  quarter  sector  computational  domain  of  an  RL-10  rocket 
nozzle  exhausting  over  a  backward  step  into  a  cylindrical 
shroud. 
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Figure  20.  Longitudinal  mesh  plane  for  the  problem  of  Figure  19. 


Figure  21.  Pressure  contour  plot  for  the  mesh  plane  of  Figure  20 


Figure  22.  Velocity  vector  plot  for  the  mesh  plane  of  Figure  20. 


X 


Figure  23.  Mach  contour  plot  for  a  cross  flow 
mesh  plane  in  the  shroud. 


computed  boundary  point  operator  differencing 
to  frozen  data  in  adjacent  patches;  (solid 
line)  effectively  continuous  grid  method. 


Figure  25.  Boundary  condition  interpolation  procedures 
poiatlon  points  _X_  ,  interpolant  data  points  Q  . 


Solution  points  _0_  ,  inter- 


25a.  Triangle  interpolation  for  overset 
mesh. 


25b.  Line  interpolation  for  composite 
mesh. 


i 


.-V'T 


25c.  Line  interpolation  for  overset 

mesh.  Straight  line  extrapolation. 


25d.  Line  interpolation  for  overset 
mesh.  Extrapolated  lines  turned 
in  composite  mesh  analog. 
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Figure  26.  Test  of  stability  of  implicit  numerical  procedure  exchanging  frozen  bound 
ary  point  data  between  patched  grids. 

26a. Two  patch  computational  mesh,  line  26b.  Density  contours  for  mesh  26a. 
interpolation. 


Figure  27.  Pressure  contours 
computed  on  101  x  101  point 
uniform  mesh. 
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figure  28a.  Shock  aligned  patched  grids 
for  the  inlet  problem. 


figure  28b.  Pressure  contours  from  first 
order  solution  on  patched  grid. 
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Figure  28c.  Pressure  contours  from  sec¬ 
ond  order  solution  on  patched  grid. 


Figure  28d.  Sketch  of  patched  adaptive 
mesh  topology  concluded  from  present 
results  to  be  effective  for  capture  of 
flow  structure  of  the  inlet  problem. 
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Figure  32b  Mach  contours  for  Fromm's  scheme 
with  adapted  Chakravarthy-Osher  limiter. 
Scheme  Ilia. 
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Figure  34.  Segmented  composite  mesh  coordinate  topology. 


Figure  35.  Mach  contours  on  a  crossflow  coordinate 
surface  intersecting  the  backward  facing  step  outside 
the  outboard  engine  nozzle. 


Figure  36a.  Mach  contour  plots  in  the  (outboard  engine)  azimuthal 
coordinate  plane  intersecting  the  LI  corner  (Figure  34). 


Figure  36b.  Mach  contour  plots  for  azimuthal  plane  intersecting 
the  small  step  at  the  upper  symmetry  plane. 


Figure  36c.  Mach  contour  plots  for  azimuthal  plane  intersecting 
the  maximum  backward  facing  step  at  the  intersection  of  the  shroud 
wall  and  the  upper  symmetry  plane. 


Figure  36d.  Mach  contour  plot  in  the  symmetry  plane  through  the 
inner  and  outer  nozzle  centers. 
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Figure  40b.  Detail  showing  hybrid  mesh  topology  at  corner  wher 
overset  wraparound  patch  provides  continuity  of  the  viscous  bo 
layer  communication  to  the  base  region. 
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Figure  41.  Mach  contour  plot  for  the  solution  on  the  hybrid  n 


O  First  order  multiple  grid 


Figure  42.  Comparison  of  experimental  measurements  and 
computed  base  pressure  d’ - 1  >■  ihut  ions  for  the  micom  tactical 
mi ss i le  model  . 


